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Abstract 

New methods for integrating systems of stiff, nonlinear, first order, ordinary differential equations 
are developed by casting the differential equations into integral form. Nonlinear recursive relations are 
obtained that allow the solution to a system of equations at time t + At to be obtained in terms of 
the solution at time t in explicit and implicit forms. Examples of the accuracy obtained with the new 
technique are given by considering systems of nonlinear, first order, differential equations which arise in 
the study of unified models of viscoplastic behavior, the spread of the AIDS virus, and predator-prey 
populations. In general, the new implicit algorithm is unconditionally stable, and has a Jacobian of 
smaller dimension than that which is required by current implicit methods, such as the Euler backward 
difference algorithm, yet it gives superior accuracy. The asymptotic explicit and implicit algorithms are 
suitable for solutions that are of the growing and decaying exponential kinds, respectively, whilst t lie 
implicit Euler-Maclaurin algorithm is superior when the solution oscillates, *.e>, when there are regions 
in which both growing and decaying exponential solutions exist. 


1 Introduction 

In many scientific fields, systems of first order ordinary differential equations are used to model physical 
processes of interest. These equations are often coupled, mathematically stiff and/or nonlinear, thereby 
requiring special numerical algorithms for their solution [1]. Solutions for such systems of equations typically 
require an excessive number of very small time steps in order to retain both stability and accuracy when 
integrated with an explicit forward difference method — like the forward Euler method, which is conditionally 
stable — even when a self-adaptive time stepping scheme is used. The implicit backward Euler method 
provides a stable integration algorithm that is asymptotically accurate for both very small and very large 
time steps; however, it loses accuracy over some intermediate range of time step size [2, 3], because the 
first order Taylor expansion of the Euler method cannot simulate the exponential solution of the first order 
differential equation. 

Alternative solution techniques are developed herein, where the differential equations are cast into inte- 
gral form. The integrands of these integrals are then expanded into Taylor or Euler-Maclaurin series, and 
integrated term by term. Expanding about the lower limit of the integral, t , leads to an explicit formulation 
whose stability for exponentially decaying solutions is unconditional in the linear approximation and in the 
exact solution, but it is only conditionally stable for quadratic and higher approximations. In order to retain 
accuracy, even in the unconditionally stable cases, the time increment must be monitored in the explicit 
formulation. The explicit algorithm provides an asymptotically accurate representation when the solution 
is of the exponentially growing type, but is not asymptotically accurate for decaying exponential solutions. 
On the other hand, expanding about the upper limit, t + At, leads to an implicit iterative formulation 
whose stability for exponentially decaying solutions is unconditional. The implicit algorithm provides an 
asymptotically accurate representation when the solution is of the exponentially decaying type, but is not 


asymptotically accurate for growing exponential solutions. The explicit and implicit solutions are therefore 
complementary in providing correct asymptotic representations for growing and decaying exponential so- 
lutions. Expanding about both limits leads to an implicit Euler-Maclaurin formulation whose stability is 
conditional, but whose accuracy is often better than either of the preceding approaches, especially when the 
solutions are oscillatory in nature and contain both exponentially growing and decaying components. 

The explicit algorithms can be resolved using a self-adaptive time stepping technique, whilst the implicit 
algorithms can be resolved using Newton-Raphson iteration, whose convergence is quadratic. 

A method similar to that which is proposed here was developed by Dennis [4] to integrate homogeneous, 
first order, ordinary differential equations with exponential type solutions. He was thus absolved from 
evaluating the integral contributions in the nonhomogeneous terms; a task which forms the bulk of the 
present paper. 

The paper begins with some discussion on notation. This is followed by an overview of the various 
algorithms derived herein for the integration of first order ordinary differential equations; this is a synopsis 
of the paper. The reader who is only interested in the application of these algorithms may skip ahead to the 
last section where four examples are addressed. Sections 4-7 give the derivations of the algorithms presented 
in the overview, section 3. Section 8 provides solution methodologies for these algorithms. The last section 
applies the results of this paper to four examples: a first order ordinary differential equation whose solution 
is known, thereby allowing the accuracy of the various algorithms to be assessed; a population model for 
the spread of the AIDS virus; the classical predator/prey population model of Lotka and Volterra; and a 
viscoplastic model used by the authors for the description of metallic deformation at elevated temperatures. 


2 Notation 


In this paper we address a general system of N t first order, ordinary differential equations of the form 


Xi T i U i = jV i (for i = 1,2,..., AT; no sum on z), 


(i) 


where the Xi[t ] are the N independent variables to be solved for. These independent variables may be scalar, 
vector, or tensor valued; for example, in our native research field of viscoplasticity, N = 13 typically, where 
two of the variables are symmetric second rank tensors and the third is a scalar [3, 5, 6, 7]. The dot £ ’ ’ 
is used to denote differentiation with respect to time, t. We choose time to be the dependent variable for 
illustrative purposes, as it so often is in physical applications; however, this is not a necessary restriction on 
the theory. The square brackets [•] are used to denote Tunction of’, and are therefore kept logically separate 
from parentheses (■) and curly brackets {•} which are used for mathematical groupings. 

The parameters iU\ [Xfc[2],t] and ,-Vi [Xk[t] t t] are, in general, functions of the variables Xk[t] and t. If 
neither parameter is a function of the independent variables, X *, then the system of equations is said to be 
linear; otherwise, it is nonlinear. To simplify the notation, we use square brackets containing time to stress 
the dependence of any variable on both AT*[f] and t. For example, we write the integral AI[At] to denote 
that the integral A7 depends on either JV * [t + Af] or and At. Similarly, iUi[t 4 Af] and %V\\t 4- At] 

denote the values of the parameters iU\ [AT[f 4 A t] } t 4 At] and iV\ [Xk[t 4 At],t 4 At] at time£4A£. These 
parameters have the following notations associated with them: 


iU[t]= f iU\ [r] dr, 

Jr-O 

iU2[t] = iU { [t] } 
iU 3 [t] = iU x [t ], 

etc., 


iV[t] = f Mr) dr, 

J 7 = 0 

Mt] = ,-viM. 


where = ^[O] = 0 and t > 0. Only for the parameters U and V are subscripts used to denote the 

order of differentiation, unless explicitly stated otherwise. 

The asymptotic, implicit, Taylor algorithms are best suited for the case where iU\[t] > 0 and t € [0,oo), 
whilst the asymptotic, implicit, Euler-Maclaurin algorithms are suitable for those cases in which the sign of 
iU\[t] is arbitrary or changes during the integration. 
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3 Overview 


This paper transforms a first order ordinary differential equation into an integral form, and then exploits 
integration algorithms of varying degrees of accuracy for effecting the quadrature. Implicit solutions are 
based upon Taylor expansions of the integrand about some future time, t -f At. Explicit solutions consider 
Taylor expansions of the integrand about the current time, f. And the implicit Euler-Maclaurin solutions 
are derived from Euler-Maclaurin expansions of the integral about both the current time, t , and the future 
time, t + At . The degree of accuracy of a particular solution depends upon how many terms are kept in the 
series expansions. 

In many of the derivations that follow, and without loss in generality, it is useful to consider a single, 
first order, ordinary differential equation, x.e. 

X + U x X = V i, (2) 

instead of the system of equations given in (1), thereby removing some excessive notation. Only in the 
last two sections of the paper is it necessary to distinguish between the various equations that make up the 
system, and hence the more detailed notation of (1) must be used. 

Integrating this first order differential equation, one acquires the well known solution 

X[t + A t] = X[t) + jT +At e -m+±t}-vm) wp (3) 

which is written here in recursive form, and which corresponds to equation (14) in the text. A recursive 
solution is particularly advantageous for numerical applications. The nonrecursive form of equation (3) 
was developed by Bernoulli [8, 9] in the seventeenth century*, and will be referred to here as the Bernoulli 
solution. 


3.1 Bernoulli Solution 

Whenever U\ and V\ are both constant valued, i.t. independent of both X and t t the differential equation (2) 
is linear and the recursive integral equation (3) can be integrated exactly resulting in 

X[t + At] = X[t] e - Ul At + (1 - e~ u ' At ) . (4) 

Ol 

This is an exact integration of a linear , first order, ordinary differential equation with constant coefficients. 
The purpose of this paper is to obtain Bernoulli type solutions (both approximate and exact) for nonlinear , 
first order, ordinary differential equations by expressing them as linear equations with variable coefficients. 

3.2 Linear Solutions 

After applying a sequence of three Taylor series expansions to the recursive solution (3), we obtain the 
following implicit, 

X[t + At] ~ X[t] + Yl j* + f 1 _ g-Uilt+At] aA (expand about t + At), (5) 

and explicit, 

X[t + At] ~ A'[t] e - y »w + Ml _ e -^i[i] 4ij (expand about t), (6) 

iinear 1 approximations for the integration of a general, first order, ordinary differential equation. Similarly, 
after applying a sequence of three Euler-Maclaurin series expansions to the recursive solution (3), we obtain 
the following ‘linear’ trapezoidal approximation, 

X[t + At] ~ x[t]e-i (t/ ‘ [tl+c/l[,+A,1) ' A ' + | ^i[t]e-> (t/,[ ‘ 1+t;,[ ‘ +A<1) ' A ‘ + Vi[< + At]) • At, (7) 

*The solution obtained by John Bernoulli in 1697 was expressed as a quadrature since the integral of dz/z in the form of a 
logarithm was not generally known until later that same year. 
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for the integration of a general, first order, ordinary differential equation, where this expansion is equally 
weighted about the two limits of integration, i.e. t and t + At. 

These three solutions correspond respectively to equations (19), (21), and (25) in the text. We refer to 
these equations as linear, even though they are nonlinear in the solution variable because the three series 
approximations used in their derivation were each truncated after their linear term. The fact that U\ and 
Vi are not constant valued, in general, is the distinguishing feature between the linear approximations (5) 
and (6) and the Bernoulli solution (4) to the first order differential equation (2). 


3*3 Quadratic Solutions 

A more accurate representation of the integrated solution to (2) is arrived at by retaining the higher order 
derivatives in the three series approximations. Upon truncating the Taylor expansions after their quadratic 
terms, we obtain the implicit, 


X[< + A<] ~ X[<]e~( l/ ‘ t<+A * 1 a< 5 ) + 


/ U 2 [t + At) y f (2>»)!K 1 [< + A<] ( _ Uilt+At] . At ^ (Ui 
^V2tM< + A<]V 1 Ui[t + At]{ 1 2^ 


[< + At] •AO" 


_ (2n + l)! Vrt+At) _ ^ \ 

n! U\[t + At] 2 y ^ it! )]' 


( 8 ) 


and explicit, 

X[( + A<] ~ X[<]e-( l/ ‘W A ‘ + i t/ >W A,, ) + 


(2n+l)! I 

n! U^tY 


f' ( Ml X / ( 2 ") ! m ( 

toKmtrJ \ n\ J mi 

- - e-^W a« j | f 


_ _-y 1 w 


E(-d 


* (W_A<) 

k\ 


( 9 ) 


‘quadratic’ approximations for the integration of a general, first order, ordinary differential equation. Like- 
wise, by truncating the three Euler-Maclaurin expansions after their quadratic terms, one obtains the fol- 
lowing ‘quadratic’ Euler-Maclaurin approximation, 


X[< + A<] ~ XW e - A ' t [ A, ] + i(K 1 [i]e- A4, f A, ] + l / i[< + A<]).A< + 

+ n (V2[<]e -A$[A<1 - ^ 2 [< + A<])-A< 2 + 

+ T2 (VxW e- A4, f A< J + Ui[t + Al]) - i(2 U 3 [t] + U 2 [t + A<]) A< - l^W'Ai 2 ) - 

- U x [t + &t]V x [t + A<]}-A< 2 , (10) 


for the integration of a general, first order, ordinary differential equation. Here 


A*[A1] = + (/,[< + A<])-A< + ±(U 2 [t] - U 2 [t + A<])Ai 2 

is the argument of the exponential used above. 

These quadratic approximations correspond to equations (37), (41), and (44) in the text. The quadratic 
implicit and explicit approximations reduce to their linear counterparts stated above whenever and 
V 2 are both zero valued. However, the quadratic Euler-Maclaurin approximation does not reduce to its 
linear counterpart under similar conditions, because of the U\ and Vi terms in the curly brackets {•} of 
equation (10). Only when U 2 and V 2 are both zero valued and A* -C 1 does the quadratic Euler-Maclaurin 
solution reduce to its linear counterpart. This is because the correct asymptotic solution of the differential 
equation (2), i.e. V\/U \ , resides solely in the linear term of the Taylor approximations, but it does not 
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reside in the linear term of the Euler-Maclaurin approximation. Higher order terms are needed in the 
Euler-Maclaurin solutions to obtain a good asymptotic solution to the differential equation. 

The infinite sums on n in (8) and (9) must be truncated in numerical applications. Their rates of 
convergence are governed by the ratio U 2 /U j*, and by the expressions in the last parenthetical bracket (■) 
appearing after the sum on n. Both of these expressions go to zero as n goes to infinity. Therein the sums on 
m represent the first n terms in an exponential series. The difference between the sum and the exponential 
goes to zero with increasing n even more rapidly than the expressions themselves. 


3*4 Exact Solutions 


Finally, exact solutions for the integration of (2) are obtained by retaining all the derivatives in the three 
Taylor series expansions that are used in its derivation. For the implicit solution we obtain 

X[< + A<] = X[<]e- A,, [ A ‘J + 

( i-)"-* ( n + 2m )- ( K»-*+i[* + At] \ 

h*hh ' k[ m!(n-i)!Ui[< + A<]" +2m+ V 


x ( 1 - e-*M*+ A ‘]- AI J2 

\ p= 0 


n+2m (Ui[t + At ]- Atf 


y 


whilst for the explicit solution we obtain 
X[t + At] = X[t] + 

+ C“ 


(A«[A(]-t/,[(] At) V' V'V'fc , / 1 \r»+2m (" + 2m ) ! f Vn-k+1 


[t] 


+2m + l 


r n+2m 


E (-!)' 

P-0 


m — 0 n = 0 k—0 


- e -^W At j , 


which correspond to equations (57) and (63), respectively, in the text. Here 

Atf[A t] = £(-l)" +1 + A t n 


n= 1 


(ii) 


( 12 ) 


and 

A9[A«] = E^A ( " 

n=l 

are the arguments of the exponential functions found in (11) and (12), respectively. Truncating the series 
in these solutions after their linear and quadratic terms leads to the linear and quadratic approximations 
listed above. Other higher order approximations, e.g. cubic and quartic, can likewise be derived from these 
exact solutions. As an aid to the reader interested in their development, Table 1 presents some of the lower 
order elements that define the matrices b m>n given in equations (52-54) and (60-62) for the implicit (11) 
and explicit (12) exact solutions, respectively. The subscript notation of the matrix b m n denotes summing 
indices; they are not differentiation indices. 

For the numerical analyst interested in evaluating a large number of derivatives, the explicit solution 
would be the appropriate ‘near’ exact solution to use, because it does not require iteration. The exact 
implicit solution is presented for those who are interested in constructing higher order implicit algorithms 
than the quadratic one which is presented herein. 
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Table 1: Elements of the matrix 6, 



b _m mu 3 

b4 ‘ °~l6’ b4,1 ~ 12 ’ K 2 


_ mu \ + 2 mm 


^4,3 = — 


9 mm + ^rnmm + 20 u 2 v 


i _ 36f/| f/ 6 + 216t/ 2 2 f/ 3 f/ 5 + 135C/ 2 2 6' 2 + 360f/ 2 f/|f/ 4 + 40f/ 3 4 
4,4 51840 


64,5 = 


_ i8t/ 2 3 t/ 7 + vi&mu^m + m uju^ + 252u 7 mm + wmmrn + i4oj/ 3 3 c/ 4 

181440 


, _ 5 mu, L _i^mm+4omm 

5,0 qo > 5,1 Qf? > s , 2 nro ' 


. 9u,u s + 60 mmm + 40 mm 

0 S3 . . _ _ , 


' 32’ a '‘ 96 ’ 1152 ’ 5 - J_ 3456 

L _ 9 mm + 72 mmm + 45 mm + i8ot/ 2 2 i/ 3 2 <7 4 + 40 mm 

5,4 20736 

, _ 21mm + 252 mmm + mmmrn + "wrnmm + 945 mmm + Mommm + ™ m 

5 ’ 5_ 435456 
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4 Recursive Integral Equation 

One can integrate the differential equation (2) over the time interval [0,tf + A*] by introducing an integrating 
factor, and thereby obtain 


X[t + At] = X 0 + 


rt+At I" /■!+£< 

/ exp - / 
j(= o L Jt=( 


dU[i 


dr 


■dr 


dm 

dZ 


dZ 


rt+&t 

- Xo e- u[(+At] + c -(t/[«+A«]-t/[f]) 

•/{= 0 


dZ 


dZ , 


(13) 


where X 0 = X[0]. The integral in this equation is the nonhomogeneous contribution to the solution of the 
differential equation. This integral equation can be cast into a recursive relation by dividing the interval of 
integration into two parts, viz . 

X[t + At] = Xo + f M ^ + f + “ e -(it[. +A «]-W[tD jf. 

Jf=0 j(z=t 

Substituting the identity 1 = exp[— U[t]] exp[t/[t]] into the first integral of this equation results in 
X\t + At] = Xoe-^+^ + e-^+^-^'D /' e -C"W-^Kl) + 

rt + bt 


+ 


c -(l/[«+A«]-t/K]) ^R] 


/«=« ae 

which simplifies to the desired recursive integral equation 


x[t + At] = x[t] e -m+±<]-m) + /“ +At e -(^+A«]-uK]) 


(14) 


whose derivation made use of (13) when integrated over the time interval [0,*j. 

The recursive integral equation (14), which is an exact solution to the first order ordinary differential 
equation (2), is practical only when a solution exists for evaluating the integral 


/•t+At 

A /[At] = / 


-(i/[t+A<]-tfKD dm d t 

e dZ 


(15) 


which appears in (14). This integral will be referred to as the nonhomogeneous integral, because it represents 
the nonhomogeneous contribution in the solution to the first order ordinary differential equation (2). The 
need to obtain solutions to integral (15) in order to develop integration methods for stiff, first order, differ- 
ential equations was pointed out by Hamming [10], but as for its solution, he left “the details to the reader” . 
In this paper several different solution strategies are presented for this integral. They differ in their approach 
(:.e. implicit and explicit Taylor, or implicit Euler-Maclaurin) and in their accuracy of approximation (i.e. 
the number of terms kept in their series expansions). 


4.1 Asymptotic Integration 

For the asymptotic Taylor algorithms developed in this paper, we assume that Ui[t ] > 0 for t £ [0,oo), so 
that U[t] = U\[t] dr is a monotonically increasing function of time. The nonhomogeneous integral (15) is 

then a Laplace integral [11, 12, 13, 14] where the integrand has its largest value at the upper limit, t + A/, 
and possesses an evanescent memory of the forcing function dV[£]/d£. This fading memory means that the 
solution will depend mainly on the recent values of the forcing function, and by concentrating the accuracy 
on the recent past we obtain accurate asymptotic representations of the solution. 

In the implicit solution, the integrand is expanded in a Taylor series about the upper limit, t + At, where 
the integrand has its largest value and contributes most to the integral. By retaining just a few terms in the 
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Taylor series expansion the integrand is accurately approximated where it is largest, and the neglect of the 
higher order terms is only felt near the lower limit, t , where the integrand contributes only a small amount 
to the integral because of its exponential decay from the upper limit. The neglect of the higher order terms 
in the Taylor series thus results in an algorithm that is asymptotically correct at the upper limit. Normally, 
when treating asymptotic expansions, the exponential decay of the integrand allows the lower limit to be 
replaced with zero or minus infinity to ease the integration. This is not done in the present case, however, 
so that by retaining the lower limit as t, we obtain a uniformly valid asymptotic algorithm in the implicit 
approximation whenever U\[t] > 0. When the infinite Taylor series is kept the exact integral is obtained. 

In the explicit solution, the integrand is expanded in a Taylor series about the lower limit. The neglect of 
the higher order terms in the Taylor series results in an integrand that becomes progressively more inaccurate 
as we approach the upper limit of the integral, where the contribution from the integrand is most important. 
The explicit approximation is not, therefore, a valid asymptotic representation of the integral when the 
Taylor series is truncated at a finite number of terms. However, the exact solution is recovered when the 
infinite Taylor series is kept in the integrand. This exact solution necessarily differs from the exact implicit 
solution. 

Whenever U\[t] < 0 and therefore U\t] becomes a decreasing function of time, the reverse situation occurs. 
In this case the asymptotic solution is now obtained by expanding the integrand about the lower limit, t f 
where the integrand contributes most to the integral. Now the implicit method, obtained by expanding 
about the upper limit, does not give a valid asymptotic representation of the integral. 

In the case where {/*[£] changes sign in the time interval from t to t -f A t, it may become necessary to 
determine the time £ E [t, t + At] where this occurs. We then have 

AI[At] - f df + / <+A< e -M«+A«J-i/[e]) df 

1 J J(=t d£ J (=( d£ 

If U i[£] > 0 for £ E \t y C] an d J7i[£] ^ 0 for £ E [C> f 4- At], an asymptotic representation of (15) can be found 
by expanding the integrands in both integrals about £ = On the other hand, if U\[Z] < 0 for £ E [t,C] and 
£A[£] > 0 for £ E [C^ + At], then an asymptotic representation of the nonhomogeneous integral (15) must 
be sought by expanding the integrand in the first integral about the lower limit £ = t , and by expanding the 
integrand in the second integral about the upper limit £ = t -f At. 

Many differential equations have regions where the derivative dV[(,]/d£ in the nonhomogeneous integral 
(15) grows more rapidly with time than the exponential decays, whilst others have regions where t/i[£] < 0 
and the integrand grows exponentially. In such cases, the asymptotic Taylor algorithms must use sepa- 
rate points about which the integrand is expanded. An alternative strategy is to evaluate the integral by 
trapezoidal or more accurate integration schemes. We may, for example, approximate the nonhomogeneous 
integral (15) with the asymptotic Euler-Maclaurin expansion. 

The Euler-Maclaurin expansion provides the value of the integral from a knowledge of its integrand, and 
the integrand’s higher order derivatives, evaluated only at the two limits, t and t + A*. This is an important 
consideration, because this type of integration method is implicit. If other quadrature methods — such as 
Gaussian integration — were used, the integrand would have to be evaluated at a number of intermediate 
points between t and t-\-At, and therefore implicit iterative procedures would be required at each quadrature 
point. 

4.1.1 Implicit Variable Change 

Implicit asymptotic solution strategies for the nonhomogeneous integral (15) with U\[t] > 0 are obtained by 
expanding its integrand into a Taylor series about the upper limit, t + At, where values for the variables 
are unknown, and therefore the solution technique must be iterative. To simplify the derivation of these 
strategies, it is useful to transform this integral into the equivalent relationship, 

AI[At] = - r e -M‘+A«]-tA[.+A»-,j) d V \t + At-_z] ^ (16) 

Jz = o dz 

by introducing the change in variable z = t + AZ — 
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4.1.2 Explicit Variable Change 

Explicit solution strategies for the nonhomogeneous integral (15), on the other hand, are obtained by ex- 
panding its integrand into a Taylor series about the lower limit, t , where values for the variables are known. 
To simplify the derivation of these strategies, it is useful to first multiply the integrand by the identity 
1 = exp [— C/[t]]exp [{/[<]], and to then transform the resulting integral into the equivalent relationship, 

AJ[Ai] = e -W*+A0-^W) e m+*]-vm dVl±z] ^ (17) 

Jz = 0 

by introducing the change in variable z — £ — t. 


5 Linear Solutions 

The simplest approximation that one can make to gain a solution for the first order ordinary differential 
equation (2) is a linear one. Implicit, explicit, and trapezoidal formulations are considered below. These 
solutions are referred to as being linear because the series approximations used in their derivations are 
truncated after their linear term. However the resulting formulae are not linear equations, because both U[t] 
and V[t] may be functions of X[t]. 


5.1 Implicit Approximation 

By expanding the arguments of the exponential and forcing functions of integral (16) into Taylor series — viz. 
U[t + At — z] = U[t + At] — U\\t + At]'Z + h.o.t ., where h.o.t. means higher order terms, and alike expansion 
for V[t 4* At — z ] — this integral can be rewritten as 

/•At 

A /[At] ~ / Vi[t + At]dz, (18) 

Jzzz 0 

which when integrated gives 

(' - ' -■) 

for an approximation to the nonhomogeneous integral (15). Substituting this result into the recursive integral 
equation (14) leads to the desired approximation [3, 5], 

X[t + At] ~ X[t] (l - e -u.[*+At]-A»J 1 ( 19 ) 

which is the linear, implicit, recursive solution to the first order differential equation (2). The derivation of 
this result also made use of the Taylor approximation U[t + At — At] = U[t\ = U[t + At] — U\[t + At] « At 4 h.o.t. 
in the exponential term found in the homogeneous solution, i.e. the non-integral term in the recursive 
solution (14). This is an implicit representation because the parameters U\ and V\ are both evaluated at a 
future time, t 4- At, and are therefore unknown a priori. 


5.2 Explicit Approximation 


Similarly, by expanding the arguments of the exponential and forcing functions of integral (17) into Taylor 
series — viz . U[t + z] = U[t ] 4* U\[t]-z 4 h.o.t., and alike expansion for V[t + z ] — this integral can be rewritten 
as 


A/[A<] ~ / e u >lt]* Vl [t]dz, (20) 

Jz= 0 

where the series expansion U[t 4 At] = U[t] 4- U\[t]'At 4- h.o.t. was applied to the exponential term in front 
of the integral in (17). Integrating this relationship gives 




Ui[t] 
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for an approximation to the nonhomogeneous integral (15). When this result is combined with the recursive 
integral equation (14), one obtains the desired approximation 

X[t + At] ~ X[<] (l _ e -^.[«] ; (21) 

which is the linear, explicit, recursive solution to the first order differential equation (2). This is an explicit 
representation because the parameters Ui and V\ are both evaluated at the current time, t , and are therefore 
known. 


5.3 Trapezoidal Approximation 

By using an Euler-Maclaurin quadrature to expand the integrand about its limits of integration [15], we 
obtain 

AI[At] ~ §(/[*] + f[t + At])- At (22) 

as a linear trapezoidal approximation to the nonhomogeneous integral (15). Here the function /[£] is defined 
as the integrand of (15), which has the form 


m 


-(U[t+A(]-l/[«]) 

5C 


V, [4] exp 



(23) 


The integral in (23) can also be represented by an Euler-Maclaurin expansion, thereby allowing the function 
/[£] to be approximated by the expression 

m a Vi KJ exp [- \ (Ui [£] + Ui (t + At]) (t + At- 0] , 

where once again the expansion is truncated after the linear term. Combining this result with equation (22) 
gives 

AI[At] ~ | (viM + prp + A *^ , A< (24) 

for an approximation to the nonhomogeneous integral (15). Finally, substituting this result into the recursive 
integral equation (14) leads to the desired approximation, 


X[t + At] ~ X[t] e -i(tM‘]+iM<+At]) A< + I g-Ku.pi+u.ft+AiD-At + y ^ t + Af X , Af> (25) 

which is the linear, trapezoidal, recursive solution to the first order differential equation (2). The derivation 
of this result took the argument of the exponential in the homogeneous contribution of (14) to be expressed 
as an integral — similar to the way in which the argument of the exponential in the function /[£] was handled 
above — which was then expanded in an Euler-Maclaurin series and truncated after the linear term. This is 
an implicit solution because both Ui[t + At] and V\ [t 4* At] appear in the solution, and therefore it must be 
solved iteratively using a technique like Newton-Raphson iteration, as their values are initially unknown. 


5.4 Discussion of Linear Solutions 


Comparing the linear, implicit (19) and explicit (21), approximate, recursive solutions for the first order 
ordinary differential equation (2) with the Bernoulli recursive solution (4), one notices that all three solutions 
have the same functional form, except that U\ and Vi are not constant valued in the implicit and explicit 
solutions of (19) and (21). The implicit and explicit solutions differ between themselves only in when the 
parameters U\ and V\ are evaluated, viz. at some future time, t -f At , versus the current time, /. Ilowever 
this does not imply that they have the same numerical characteristics, for they do not. 

For large time steps, At 1, with U\ > 0, the exponential term becomes small compared with 1, and 
the asymptotic expansions 


lim X[*-|-Af]x 

At— ►large 


V\ [< + At] 
CM« + A<]’ 


(26) 
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( 27 ) 




and 


lim X[t -f Af] x iVi [f 4- Af]-Af 

A/— ►large 


(28) 


are obtained for the linear implicit (19), explicit (21), and trapezoidal (25) solutions, respectively. The 
asymptotic expansion in (26) is the correct asymptotic solution to the first order differential equation (2); 
the other two asymptotic expansions, (27) and (28), are not correct asymptotic solutions. 

At the other extreme, for small time steps, At <C 1, the exponential term has a value of about unity, and 
(19), (21), and (25) reduce to 


lim X[t 4* At] = X[t] 4- Vi [t + At]’ At + h.oA ., (29) 

Af — ►small 

lim X[t + Af] = X[f] + Vi[i]-A< + h.o.t., (30) 

Af — ►small 

and 

lim X[< + A<] = X[<] + i(Ki[<] + V , i[t + A<])-A< + A.o.<., (31) 

A*— ►small z 

which correspond to the classical Euler backward, Euler forward, and trapezoidal difference methods, re- 
spectively, whenever U\X \s small compared with Vi. 

Although the difference between the implicit and explicit approximations in (29) and (30) is negligible 
for small Af, the difference between (26) and (27) can be enormous for large Af. If we start at f = 0 and 
apply a large time increment Af , the solution Vi[At]/U\[At] for the implicit approximation is asymptotically 
correct, whilst that of the explicit solution, Vi[0]/(7i[0], corresponds to a ratio of the slopes of V and U at the 
initial time, f = 0. In the case of viscoplasticity, these situations correspond to a correct viscoplastic solution 
in the implicit approximation, and an incorrect elastic solution in the explicit approximation, respectively. 
In subsequent time steps the explicit approximation will oscillate around the true solution, but it will not 
become unstable. These oscillations can be mitigated only by choosing smaller time steps. The implicit 
approximation, on the other hand, is accurate for all sized time steps for exponentially decaying solutions. 

In the case where V\ is slowly varying and higher order derivatives can be neglected, the asymptotic limit 
is V\ [f 4* At]/U\[t 4- Af], and not Vi[f 4- Af]Af/2. However, the current trapezoidal approximation provides 
accurate values of A/[Af] for intermediate time steps, and this is an important consideration in those cases 
where the derivative V\ in the integral swamps the exponential decay, or when U\[t] becomes negative and 
asymptotic methods require expansions about appropriate limits. 


6 Quadratic Solutions 

Quadratic solutions for the first order ordinary differential equation (2) are derived in this section by using 

y = U\[t + Af ] • z (32) 

as a change in variable for the implicit solution, and by using 

y = U\[t]-z (33) 

as a change in variable for the explicit solution. These quadratic approximations are derived in a similar 
manner to the linear approximations determined above, but now the three Taylor series expansions are 
truncated after their quadratic terms. Likewise, a quadratic Euler-Maclaurin approximation is derived by 
expanding the integral solution (14) with three Euler-Maclaurin series, which are then truncated after their 
quadratic terms. 
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6.1 Implicit Approximation 

If the Taylor expansions used in the derivation of the linear implicit solution are truncated after their 
quadratic terms rather than the linear ones, then (18) becomes 

A/[A<] ~ I*' e -(tM‘+A<] * a ) (^[f + At] - V 2 [t + A t] z)dz, (34) 

Jz = 0 


which is a more accurate approximation to integral (16). 

By introducing the change in variable (32) into (34), one obtains 


r u i! 

AI[At]~ / 

Jy= 0 


U 1 [t + At]At , 

e -y e \u 7 [ t+A,) z (Vi[< + A< j _ + zl dy> 

dy 


which, upon expanding the second exponential in the integrand into a power series, interchanging the order 
of summation and integration, and differentiating the change in variable (32), enables integral (34) to be 
written as 


AI[At] 


1 ( U2\t±M V 

4ln! \2{/i[f + A<] 2 


n=0 


v)[< + a<] 


U\ \t *4* At] Jy—Q 


1 


e y y 2n dy 


Noting that [16] 


f e y y n dy = n! j 1 — e 
JyzzQ \ 


- ■*> 

■i S)- ~ 


which is the difference between the gamma and incomplete gamma functions, enables the two integrals in 
(35) to be integrated, giving 


A/fAtl ~ V 

1 J " \2{/i[t + At] 2 


/JM<+A<]_y / (2n)! Vi[f + Af] 


2n 


1 _ _-lM‘+A<] A« Y' 

Vfi + At] I 


(iMi + 


ml 


(2n -f 1)! V 2 \t±M 


1 


n! U { [t + At] 2 \ 

This leads to the desired approximation 

X[* + Af] ~ X[<] e -(tM(+A<] A<-lCMt + Ai] At 3 ) + 


£/,[t+At)At yt 1 (Jh[t + A<] A<)* 


k=0 


k\ 


)}■ 


+ 


Y / U 2 [t + Af] \ n r (2n)! Vi[t -f At] / _ 
\2t/i[< + At] 2 ) | n! C/j[< + A<]( 

_ (^ n + j2[ Y^L ( 1 _ „-f/i[f+At] At Y' 

n! Urlt + At^y 


-^[t+Atl-At v (^i[< + A<] A<r 


m=0 


ml 


(Uxlt + AtyAt) 

k\ 


))• 


(37) 


which is a quadratic, implicit, recursive solution to the first order differential equation (2). This is an implicit 
representation because the four parameters U \ , t/ 2 , Vu and V 2 are all evaluated at some future time, t -f At } 
and therefore the solution algorithm must be iterative, as their values are not known in advance. 


6.2 Explicit Approximation 

By truncating the Taylor series expansions used in the derivation of the linear explicit approximation after 
their quadratic terms rather than the linear ones, (20) becomes 

A/[Al] ~ e -(ch[<] A<+*iMt]'At 3 ) j*' e (tMi] z+JtMO *’) (y,[<] + y 2 [<] . z ) dz, (38) 

Jz=0 
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which is a more accurate approximation to integral (17). 

Introducing the change in variable (33) into (38) gives 

fUilt] At 

AI[At] ~ e-(^.W-A*+iw.[i]-Ai») f e y (V^[<] + U 2 [t]-*) ^ dy, 

Jy- 0 “ 2 / 

which, upon expansion of the second exponential in the integrand into a power series, interchanging the 
order of summation and integration, and differentiating the change in variable (33), enables integral (38) to 
be written as 


[U >l 
Jy = 0 


ViM r" 

Ui[t] J y= o 


C/,[<] A< 


e v y 2n dy + 


m 

Ui[t ? 




e y y 2 " +1 dy 


Using the fact that [16] 




(39) 


(40) 


[t] A<) m _ u t [t) At I _ 


allows the two integrals in (39) to be evaluated, resulting in 

a/[a<] ~ y ( _MLV / MZlI ( y^iyntPi 

1 j " 1 »i } 

(2n + 1)! Vb[<] f 2 sr^ ( .Nib (^iM -At)* -t/.fii-At) 1 

W\to U )}' 

This leads to the desired approximation 

_X[I+AI] ~ x[(Je" (t '> wl, ' t i l,J 'l 4l ’> + 

. '-mw v ( X ! < 2 "> ! ( f, I'm tWAg _ 


At 


_ (2n + 1)! Vn[t) r^ ( 

nl 1 k\ 


,A< ) t _ P -Vx [<] 


*))■ 


(41) 


which is a quadratic, explicit, recursive solution to the first order differential equation (2). This is an explicit 
representation because the four parameters U\, V\, and V 2 are all evaluated at the current time, t, and 

therefore their values are known. 


6*3 Euler-Maclaurin Approximation 

Let us express the nonhomogeneous integral (15) in terms of an asymptotic Euler-Maclaurin expansion [15] 
in the form 

7[A<] = i(/[t + At] + f[t])-At + ( ~/ 9 )n \f ~ (/a»-i I* + A<] - h n -i[t])-At 2n + h.o.t. 

( 2n ) ! 

= §(/[« + At) + /[<])• A* - £ (/1 [t + At] - fi[i\) At 2 + £- 0 (Mi + A<] - h[t])-At A + h.o.t., (42) 

in which subscripts on /, like U and V, denote the order of differentiation. The Bernoulli numbers, B ny are 
defined as 

- ^1 = 6» -^2 = 30 > ^3 = 42 > = 30 ’ = 66> ^6 — 2730 5 e tc., 
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where these subscripts do not denote differentiation. The series is asymptotic, or semi convergent, and must 
be truncated at k terms before the series begins to diverge. The function /[£] is the integrand of equation 
(15), as defined by (23). The integral in the definition of /[£] can also be represented as an asymptotic 
Euler-Maclaurin expansion, and consequently the function /[£] assumes the form 

m = Vi fcj exp [- (| (tf, KJ + lh[t + At ]) (i + At-0 + 

+ E (^n[< + At] - c/ 2 „k]) (t+At- o 2n + h.o.t. )] 

r»=l ' ' 

= ViK] exp [- (i (Udt] + lh[t + At]) (t + At - 0 - 

- £ (U 2 [t + At] - tf 2 K])(< + At - 0 2 + m ( u *[t + At] - t^K])(t + At - 0 4 + h.o.t. )] . 

By retaining just the first terms in these expansions, we approximate the integral with a trapezoidal formula 
which involves the evaluation of the integrand /[£] at times t and t + At. Because U\\t + At] is unknown, 
the resulting formula for X[t + At] is implicit and must be resolved with Newton- Raphson iteration. 

If we retain first order trapezoidal terms, the nonhomogeneous integral (15) is approximated by (24), 
whilst the retention of quadratic terms yields the improved approximation, 

| (viW + V x [t + At ]) At + 

A (V 2 [t] e -**l*<] _ V 2 [t + At]) At 2 + 

h {ViW (i(l/i[t] + Ui[t + At]) - 1(2 U 2 [t] + U 2 [t + At]). At - ^%[t]-At 2 ) - 

— U\\t + At] Vj[t 4- At]) • At 2 , (43) 

where the argument of the exponential is defined by 

A* [At] = \ (Ui[t] + t/i[t + At])- At + i _ U2 [ t + At])- At 2 . 

Substituting equation (43) into the recursive integral (14) gives the desired approximation 

X[t + At] ~ X[t] e - A *[ A ‘] + i (viW e -A ’ I '[ A( ] + V)[< + At]) At + 

+ 12 (v 2 [t] e -A *^ A< J — V 2 [t + At]) At 2 + 

+ n{Vi[<]e" A * [A() [<] + Ux{t + At]) - i(2(/ 2 [t] + U 2 [t + At]) -At - ^t/ 3 [t] At 2 ) - 

- Ui[t + At]Vi[f + At]} At 2 , (44) 

which is the quadratic, Euler-Maclaurin, recursive solution to the first order differential equation (2). This 
quadratic solution has an advantage over the implicit and explicit quadratic solutions, in that it does not 
require the retention of an infinite series. In addition, it contains information about the parameters U \ , 

Vj, and V 2 at both ends of the interval, t and t T At; whereas the implicit and explicit Taylor solutions 
contain information at either one end or the other, but not both. 


AI[At] ~ 
4* 
+ 


6,4 Discussion of Quadratic Solutions 

For large time steps, A* > 1, with \Ui[t + At}' At 2 < U\[t -f At] -At, the exponential term becomes small 
compared with 1, and the asymptotic expansion 


lim X[t-bAt] 

At— ►large 


( ^ 2 [t 4- At] \ / (2n)l V\ [t -f At] (2n 4- 1)! V^\t + At] \ 
^ V2f/i[t-f At]V \~^T Ui[t + At] n\ U x [t 4* At] 2 ) 
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V\ [/ 4* A/] / V\ [i + At] U 2 4~ At] — U\ [t + At] V 2 [t + At] 

Vi[t + At] + V C^i[t + At] 3 


x 


^(2n-l) 


n = l 


(2n — 2)! 

(»“!)! 


/ t/ 2 [t + At] V 
V 2 C/ x [< + At] 2 / 


x 


(45) 


is obtained for (37). As the solution goes to saturation — a property of first order differential equations with 
exponential type solutions where U\ > 0 — this asymptotic expansion reduces to that of the linear case, (26), 
because the higher order derivatives f/ 2 and V 2 vanish there. Therefore, the added benefit of retaining the 
higher order derivatives U 2 and V 2 manifests itself in improved accuracy, but only in transient domains; the 
linear implicit solution is of equal accuracy in the neighborhood of saturation. Like the linear case, physically 
meaningful asymptotic expansions for the quadratic explicit (41) and Euler-Maclaurin (44) solutions do not 
exist when the arguments of the exponential functions are positive. 

For small time steps, A< < 1, the exponential term has a value of about unity, and (37) reduces to 

Jau*i ,+A< i = + A( ) - griSti}} <«> 

for the implicit case, whilst (41) reduces to 

A< 1 = + (*M - t|) A ' + t " <"> 

for the explicit case. These limits add the quantity -V^-AZ/t/i to the analogous limits given in (29) and 
(30) for the linear approximation. In contrast, the limit in the small for the quadratic Euler-Maclaurin case 
is identical to the linear case given in equation (31). 


7 Exact Solutions 

Implicit and explicit, exact, recursive solutions for the first order ordinary differential equation (2) are derived 
in this section. From these solutions, one can construct cubic and higher order, implicit and explicit, approxi- 
mate solutions. A similar procedure for constructing cubic and higher order Euler-Maclaurin approximations 
follows along the path of development given earlier in the quadratic Euler-Maclaurin solution. 


7.1 Implicit Solution 

From an expansion of the integrand into an infinite Taylor series about the upper limit, the integral in (16) 
may be written in the form 


AI[At] = [ A * exp \ — U\[t + At]- z\ exp f](-l 

Jz = 0 1 J [n = 0 \n-\-Z)\ 


> 4- At] 


X ^(-1)P '^ L V z p 


< P-0 


P ! 


j dz. 


For brevity we set 


b — (_n n ^2+ilLi^l anc j c — ( op ^ 

n-( ^ (n + 2)! c p-( 1) p , 


(n + 2)! 

and expand the last exponential in (48) as a power series, so that 


r&t if 00 r 2m / 00 \ m °° ) 

AI[At] = j^exp[-U l [t + At]-z} j 


( 48 ) 


(49) 


( 50 ) 
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The first task is to evaluate the central sum in (50). If we put 


= E 6 " 2 "’ 

n=0 

then the Cauchy product g[z]g[z ] is given by 

CO 

g[z] 2 — (bob, + bibv-i + • ■ • + 4- 6|/6o) z v » 


v -0 


which is obtained by collecting the coefficients of z v in the two product series. Following Knopp [17], we 
call the coefficient of z v in the preceding series b 2)V) so that 

oc 

a w 2 = E^. 


u = 0 


where 


Then, proceeding in sequence, 


b2 t u — bob , -j- + • • • + b v -\b\ -f b v bo. 


where 

and in general, 


AW 3 = E 63 ^"’ 

y = o 

&3,i/ = bob 2t v 4* bib 2t ,-l + h — 1^2,1 + &i/&2,0, 


i /=0 

where the 6 mji/ are obtained from Glaisher’s relation [18], 

bm,v — ^O^m— 1,1/ 4" b m _ i ^ _ i 4- 1- b, — i b m _ i i -|- b,b m _i o 

v 

— Y J 


(51) 


* = 0 


Substituting into (51) from (49) we see that the are obtained recursively from the lower order b m _ Y , 
according to the relation 

A __ V^/ t xJfc ^Jk+aP 4- Af] 

bm,n - _ rTZTv' 6m ”bn- 

*=0 


Because 


and 


(*4-2)! 

0 

^[^]° = EZ = i 

|/=0 

oo oo 

ffW 1 = E 61 ."^ = ’Y J b v z v , 


-*• 


(52) 


i/=0 


i /=0 


we must require that 
and that 


6 0( o = 1 and 6 0 ,n = 0 (for n = 1, 2, . . .), 

6,, n = b n = (fern = 0,1,2,...). 


(n + 2)! 

Tabulated values for the coefficients 6 m n up to m = 5 and n = 5 can be found in Table 1 . 


(53) 

(54) 
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The expression within the braces in the integrand of (50) may now be written as B[z \ , where 



By using the same techniques we may write the Cauchy product 


OO 

oo 

/Ev 1 ’ 

oo 

— ^ (^m p 0^n d" b m \C n —i + • ■ ’ + b m) Yi — \C\ d” % 

u—0 

‘a 

II 

o 

n= 0 



oo 

- d z n 

— / u m,rj * t 
n=0 


which when combined with (49) results in 

, _Y^. - - - YY iyi-Jt , 

Um,n — / j bm,k Cn — k — / t) / £\j Vm,k' 


k = 0 


k = 0 


(55) 


The braces may now be written as 


oo 2m °° 00 w i 

- E - E £ *” +2 ” 


„ m. 
m = 0 n=0 


m=0 n = 0 


Introducing the change of variable (32) into (50), along with the expression for B[z] ) and interchanging 
the order of summation and integration, allows the exact nonhomogeneous integral (15) to be written as 


/•LMt+Afl At co oo , / y 

A,lM] = L 


n + 2m 


dy 


m=0 n = 0 


Ui[t + A<]' 


(56) 


Recalling the integral solution (36), and substituting for d m n from (55), enables this integral to be integrated 
yielding 

ATTA/i _ V V V /) < i\ n-t (» + 2m)! / V^-t+i^ + Af] N 

A/[A<] _ !) m!(n — Jfc)! + Ai]" +2m+1 / 

m=0n = 0fc=0 V ' X 1 


n sr^ ( u i[* + AQ-A*) P 


x 1 - 


£ 

P= 0 


p] 


in which the coefficients b m * are obtained recursively from equations (52), (53), and (54). 

An exact solution for the differential equation in (2), as given in (3), is therefore provided by the relation 

X[< + A<] = + 

+ 


Y'Y'Y^/) ( iy»~* ( n + 2m)! f Vn-t+i[< + At] \ 

l > ra!(„_t)! \U 1 [t + At] n + 2m + 1 J 

x A _ e -c/,[(+At] At + At] AQ p \ ^ 

V p=o p ‘ ) 


(57) 


which is an implicit recursive solution, in which the coefficients b mt k in Table 1 depend on the derivatives 
U p [t H- A t] through equations (52-54), and where 


Ai?[A<] = ^(-l) n + l t/n ^ A ^ A<" 


nzzl 
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defines the argument of the exponential in the homogeneous contribution to the solution. This solution is 
unconditionally stable provided that U[t] is a monotonically increasing function, re. At?[Af] > 0. If on the 
other hand U[t] is a decreasing function of time, then the solution grows exponentially and the size of the 
time step must be monitored. For this case, as stated in section 4.1, it is the explicit solution, developed in 
the next subsection, which provides the asymptotically accurate expansion. 


7,2 Explicit Solution 

By expanding the integrand into an infinite Taylor series about the lower limit, the integral in (17) may be 
written in the form 


A/IA,] = .-M £ „p [c, H .,] exp [g *”«] ( £ 5±}H 


z p dz , 


where 


A»[Ai] = £ 

n = l 


Proceeding as in the previous section, we put 

u n+2 [t] 


b n = 


(n + 2)! 


and c„ = 




for brevity, and by the change of variable given in (33), obtain 


r Ui[t]-At oo oo , / \n+2m 

AI[At] = e~ A *M / V; (-^ ) 

4=o m! 


dy 

Ux[t]' 


in which 


with 


j _ f' Kt-H-iW i 
dm -" - 2 ^ („ _ k y_ 


*=0 


b 0 0 — 1 and t 0 ,n = 0 (for n = 1,2,.. .), 

U n+ 2[t] 


bl n — b n — 


(for n = 0,1,2,...), 


and 


(n + 2)! 
h -\^hh _v 

— / j vJfcPm — l,n — k — / 4 ^ 2^| l t n — fc* 


Tabulated values for the coefficients b m>n up to m = 5 and n = 5 can be found in Table 1. 
On effecting the integration via (40), we obtain 


(58) 


(59) 


(60) 

(61) 


(62) 


A/[A<] = e-^^- 

/ n+2m 


E <-» 

P = 0 


C/,[l] A«) V'f'rj / ,\n+2m (» + 2m)! / + 

hhh { m! (»-*)! UiW B+2m+ V 


where the coefficients 6 m fc are defined in equations (60), (61), and (62) in terms of the derivatives and 

where A0[A<] is given by (59). 
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An exact explicit solution to the differential equation (2) is therefore given by 
X[< + A<] = X[<]e- A ^ A( ] + 

, -(A#[A<]— tM<] A,) y* y' y' b / ^n+am ( n + 2m )! f ) 

+ e l > m \(n-k)\ Vt/![<]"+ 2m +V 

m=0n=OJt = 0 V y X lJ 


f n+2m 

E 

S P—0 


x l E ( p-IAtM-Al 


pi 


( 63 ) 


where all the derivatives, U p [t] and V p [t], are evaluated at the current time, t , and tabulated values of b m> * 
are given in Table 1. Since A#[A2] represents an exact infinite series, this solution is unconditionally stable 
provided that U[t] is a monotonically increasing function of time, i.e. A0[At] > 0. When truncated to a 
finite number of terms, the solution is conditionally stable and the size of the time step must be monitored 
for accuracy and stability. 


8 Solution Methods 

The previous sections of this paper contain derivations for implicit, explicit, and Euler-Maclaurin solutions 
with varying degrees of accuracy for the first order ordinary differential equation (2), as described in (3). 
This section presents algorithms for solving a system of first order ordinary differential equations (1), and 
we illustrate these procedures by using some of the integration methods derived above. 


8,1 Implicit Algorithms 

A technique like Newton-Raphson iteration must be used whenever a system of N , first order, ordinary 
differential equations is to be solved simultaneously using an implicit integration algorithm. For illustrative 
purposes, consider the following solutions for a system of N equations, i.e. 


X i[t + At] ~ Xi[t] + 


iVi[t + At] 
i U\ [tf + At] 


^ _ e -,tMt + Ai] At 


). 


(64) 


and 

Xi[t + At] ~ Xi[t] e-*(^‘W+^*t‘+ A *]) A* + i ( 4 Vi[<] A( + ^[t + At]) At, (65) 

for i = 1,2 These are the linear, asymptotic (19) and trapezoidal (25), implicit solutions for the 
recursive integral equation (14) when written for a system of equations. A similar solution strategy to that 
which follows can also be constructed for the higher order implicit solutions. 

In many cases, not all of the functions iU\[t + A<] are independent. For example, in the case of uni- 
fied viscoplastic constitutive equations, which are usually comprised of thirteen equations, the number of 
independent functions is usually two or three (see Example 4). Let M denote the number of independent 
functions comprising the set iUi[t-\- At] for i = 1 to AT, with M (E [1, JV]. Also, let Q a [t + A*], for a = 1 to M , 
denote either the independent functions iU\[t + At] in (64) or the independent functions ^(*{/i[f] + %Ui[t + Ai]) 
in (65). The solution vector X and the forcing function vector V\ can therefore be written in forms such as, 


{X 1 ,X 2 ,...,X i) ... 1 X / v) T = {Xi,X 2) . 

• • » X a , . . . , 

Xa/} T 


/ y 1 V 2 y^i Y 1 y 2 Y N * 

1 j • * 

v y v V 

y 1 y2 
. , A 0 , A a , . 
s. - 

y Nc 

V 1 V 2 vNm 

• > > A Af> * ■ • » ^ M 

associated associated 

with Q\ with @2 

V 

associated 
with Q a 

associated 
with qm 


where the vectors X Q and a V\ comprise the set of all components X{ and \> respectively, over the range of 
i = 1 to N, that are associated with the variable Q a . Each of the M component vectors, X a and a Vi , are 
comprised of N Q elements, where = X. 
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We may now write (64) as 


and (65) as 




XM = X a [t}e~ e °* t + \ 


A \x y [e„] 1 


3 

s 

1 

1 

c£ 

(66) 


(67) 


where the dependence of X a [t + At] on g v [t + At] is explicitly stated, and the dependence on t + At is 
implicitly assumed. In general, the g a are functions of the X y [t + At] which are in turn, via (66) or (67), 
functions of the g u [t 4* At], and therefore 


Qa 


^7 [t?*]] “ Qa [Qv] 


(for a, 7, v — 1 to Af), 


or equivalently, 


i a \Qv\ — 0) 


(68) 


where 

■F a [Qu] = Qa [Qu] ~~ Qa- 

If {&»}* i s fbe ^ th guess for vector g a , then the true solution satisfying (66) or (67) may be written as 


Qa — H" C a} 


where the correction vector, c a , is the amount by which the true vector differs from the guessed value. 
Inserting this definition into (68), and expanding the resultant by Taylor’s theorem while retaining only the 
first order terms leads to 


a — •F a }x 4“ Cy] 

= 0. 


A solution to the resulting linear system of equations, 


M 


T Jap Cp = -T a [{ev) a] (for a = 1 to M), 
P=1 


for the unknown correction vector, c^, with Jacobian, 


J a0 — 




0/5 d{ 6 p} x ’ 


(69) 


(70) 


can be obtained by Cramer’s rule (Af < 3) or Gaussian elimination (Af > 4). This leads to an improved 
value for the solution vector, {£ a }*_|_i, where 


Kh+1 — {#<*}* d" C OT1 


(71) 


which is iterated until the contribution from c a becomes negligible. In order to retain algorithmic stability, 
it is useful to bound the size of the correction vector so that 

iicii<THf>m 


(we typically set T = 0.25) where 

IMI= V / c 1 2 + c 2 2 + --- + c A / 2 and ||e|U= + {02}* + • • •+ Xqm)\ • 
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That is, if 

||c||>T||ff|| set Vae{l,M} 

in equation (71). This slows down the rate of convergence, but it keeps the procedure stable and therefore 
free of oscillations. 

If possible, the Jacobian should be determined analytically; if not, it can be acquired numerically, but 
often at a greater numerical expense. Because Newton-Raphson iteration has quadratic convergence, the 
number of iterations, A, required for convergence is usually few in number. One reasonable initial guess 
to advance the iteration process to the next step is to use the values from the last time step, i.e. set 
£ar[^ “1" A^]a = 1 = £cr[^]- 

The advantage in using this integration approach over classical backward difference methods — like the 
backward Euler method — when the system of equations happens to be a system of tensor valued equations 
with M N is evident, because the Jacobian is now of size M x M instead of N X N . Classical approaches, 
on the other hand, iterate one scalar valued parameter for each spatial dimension, resulting in a much larger 
Jacobian of size N x N, In the authors’ native research field of viscoplasticity, there are usually three, 
simultaneous, first order, ordinary differential equations to be solved for; two are symmetric, second rank, 
tensor equations and the third is a scalar equation [7, 19]. This results in a spatial dimension of thirteen, 
yet only three scalar valued parameters need to be determined by iteration in our approach, as opposed to 
thirteen in the usual backward difference and other implicit approaches. 

This method has been used by the authors [3, 5, 6] in the solution of several different inelastic material 
models with great success. It has also been implemented into nonlinear finite element packages [6, 20] where 
preliminary results from inelastic structural analyses show the method to be stable, equivalent or superior 
in accuracy, and requiring about the same to 1/5 the cpu time when compared with the well known solution 
techniques of: Euler forward difference method using a self-adaptive time step procedure [6]; second and 
fourth order Runge-Kutta methods using self-adaptive time step procedures; and the trapezoidal backward 
difference method using Newton-Raphson iteration [20]. Because these were feasibility studies, their primary 
emphasis was not on writing efficient code, but rather to determine the feasibility of using the asymptotic 
integration method in finite element applications. Further improvements in computation times are expected 
as we become more knowledgeable about writing efficient code for the algorithms. 

8,2 Explicit Algorithms 

Different needs demand different solution techniques, and Newton-Raphson iteration will not always be the 
technique of choice. This technique is ideally suited for large computer codes where large time steps are 
needed to minimize the computation time. But it is not the best choice when small time steps are required 
in order for one to gain a detailed picture of the response. For applications of this type, explicit solution 
techniques work best. 

The simplest integration algorithm considered in this paper is the linear explicit method given by 

Xi[t + At] = Xi[t] e-' u 'W ** + (l - e-' u 'W A ‘) (for » = 1,2 (72) 

In order for this approach to give accurate and stable results, it is necessary to keep the time step small. 

A more accurate integration algorithm than (72) given above is the quadratic explicit method defined by 

Xi[t + At] = 

+ 


for i — 1,2, . . N. The infinite sum on n found in (41) is truncated here to q terms, which must be set by 
the user. The time step must also be kept small for this method, too. 


Xi[t]e-^ Ul[i] ** ) + 

V' ( V f ( 2w )~ ( VVjy. 

iiUwJ 1 »! i<h[t){ko } 

} *1 )} 


(iU.it]- Aty 




(73) 
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8.3 Predictor-Correctors 


Predictor-corrector methods use an explicit algorithm to first look ahead and predict the value of the function 
at a future point. Knowing something about the future response of this function, one then uses an implicit 
algorithm to go back and make a second, more accurate, correction of what the future value of this function 
is. These methods have the advantage of being able to use the difference between predicted, X-[t + Af], and 
corrected, X{\t + At], values to establish an error, e.g. 


e — max 


Xj[t + At] - XHt + At] 


Xi[t + At] - Xi[t] 


£ — max 

i 


X<[t + At] - X'[t + At] 


Xi[t + At] 


)■ 


which one can use to control the time step size. The first error listed is useful in transient domains where 
Xi changes rapidly; it is the error accumulated over the time step, At. The second error listed, i.e. the 
relative error, is the better choice when X{ changes slowly, because the first error term has a denominator 
that becomes small under these conditions. 

A linear predictor-corrector can be constructed by considering 


Xj[t + At] = X,[t] e-‘ u 'W At + f 1 - A ‘) 

{ U i pj \ / 


(for i = 1,2, ...,N) 


for the predictor, and 

Xilt + At] = x,[t] c-- y it‘ +A ‘] (i _ c -.</|[« + A0 A«) (for i 

for the corrector, where 


1,2 


(74) 


(75) 


iU[[t + At] — iU\ 


Xj [/ -f- At] y t -f- At 


and iV([t + At] = ,V i [xj [t + At] , t + At] . 


For this linear integration method both the predictor and corrector are unconditionally stable, but the time 
step size must still be monitored for accuracy to avoid an oscillatory response. 

A more accurate predictor-corrector than (74) and (75) given above is the quadratic predictor-corrector, 
where 


M-AQ" 


.- ( lt,[t]-AI _ 


Xj[t + At] = X,[t]e- (lfA[t] i,+ 5'W i|, ) + 

+ e -i«y»w-A« a y' f »EM < ] V f ( 2n ) ! ( y*, 

_ (2n + 1)! /yp -.t/.m.Ai) 1 

is the predictor, and 

X[f + At] = X [ f ] e -(,f/;t<+A(] + Al a ) + 

+ X' ( itffi + A*] y / (2 n)! <V/[f + A/] f -MU+fitlAt y> G^ft* + At]-At) m \ 
£ 0 \2iUi[t + At]*J \ n! iUl[t + At]\ ^ ml j 

(2n + 1)! iV’[t + At] (, ,_ iJ7;r(+All . Af (, V[[t + At]-At) k \\ 

n! iUllt + At^y k\ )}' 

is the corrector, for i — 1,2, . . N. Here 

iU[ [t + At] = i U x [xj [t + At ] , t + At] and ,-Vflt + At] = ,-V, [xj [t + At] , t + At] , 


(76) 


(77) 


with 


+ At] = iU 2 Xj [t + At], t + At and ,V 2 '[t + At] = iV 2 [xj[t -f At], t + At 


The size of the time step for this quadratic predictor-corrector must be monitored for both stability and 
accuracy so as to avoid an oscillatory response. 
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9 Numerical Examples 

In this section we consider four initial value problems to explore the capabilities and deficiencies of the various 
algorithms derived in this paper for the integration of systems of nonhomogeneous, nonlinear, first order, 
differential equations. The first example is simple enough to enable an exact solution to be derived, but yet 
sufficiently difficult to tax the integration algorithms. This example is used to compare the accuracies and 
rates of convergence of the various implicit integration algorithms developed herein. The classical method of 
backward Euler integration is used as the meter stick for comparison. The second example is a population 
model developed for predicting the growth of the AIDS virus in the human population. Like the first example, 
this example possesses an asymptotic solution at large time; however, unlike the first example, its path to 
the asymptotic solution is not monotonic. This makes the numerical integration much more difficult. The 
third example is the classical Lotka-Vol terra predator-prey model. This is also a population model, but with 
an undamped oscillating solution. As a consequence, damping present in an integration algorithm — either 
under- or overdamping — has the potential of adversely influencing the predicted result. The final example 
is a viscoplastic model used to describe the deformation behavior of metals at elevated temperatures. This 
model possesses asymptotic solutions; therefore, asymptotic integration algorithms work best. Unlike the 
prior examples, its Jacobian cannot be reduced to a 1 x 1 matrix; although its dimension has been reduced 
from the spatial dimension of 13 x 13 to a 2 x 2 for the Jacobian. 

9.1 Example 1 

The nonlinear differential equation 


y + y 3 = 1 with y(0) = 0 (78) 

mentioned by Krieg [21] may be solved by the method of separation of variables, and the solution is 


t = — = tan~ 1 

V3 


1 + 2 y 


L v/3 J 


- 3 ln[l ~y] + h ln[l + y + J/ 2 ] - 


6X/3' 


(79) 


The variation of y with t is shown in Fig. 1. 

The differential equation (78) may be written in the notation of this paper as 

y + U\ • y = Vi where U\ — y 2 and V\ — 1. (80) 

Higher order derivatives of U are easily calculated from (78) as 

U 2 = lh = 2yy = 2y(l - y 3 ), (81) 

U 3 = U 2 = 2y(l - Ay 3 ) = 2(1 - y 3 )(l - 4y 3 ), (82) 


v 2 = V 3 = ■ ■ ■ = 0 . 


(83) 


The solution, according to either equation (37) or (57), up to terms including U 3 , is given by the quadratic, 
implicit, recursive relation 


y[< + At] = y[t]e-**^l + (l - + 

+ Vl[t { l _ e ~ u ' [‘+*‘1 At (l + Ui [t + At] ■ At + \ (Ui [t + A<] • At) 2 ) } + 

+ M! { i _ e -cM<+At] At ^ + Ul [ t + At]. At + i {U x [t + At]- At) 2 + 

+ i (Ui[t + At] -At) 3 + ± (Ui[t + At]- At) 4 ) ) , 
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where 

Ai?[A<] = Ui[t + At] -At - ±U 2 [t + At] At 2 . 

Let 

ei=U 1 [t + At] = y[t + At] 2 , (84) 

so that from (81) 


U 2 [t + At] = Ui [t + At] = 2 y[t + At] (l - y[t + At] 3 ) = 2 e \' 2 (l - qT) ■ 

The argument in the exponential of the homogeneous term can now be expressed in terms of e\ , 

Ai?[A<] = U l [t + At}- At- \U 7 [t + At] At 2 
= ^-At - e\ >2 (l - e? /2 )-At 2 , 

and the quadratic, implicit, recursive solution is obtained in terms of q\ in the form 

y[t + At] = i/[t]e-( ei A ‘ 3 ) + l(l_ e -^ A, ) + 

ei 

2e\ /2 (i-e 3 l /2 ) ( 

+ ^5 1 1 - e ~ ei At (i + Qi -At + | (ei -At) 2 ) } + 

+ — - 2g ‘ ( ;; fi ^ 

+ | (ei At) 3 + ^ (ei At) 4 j | . (85) 

Similarly, the quadratic, Euler-Maclaurin, recursive solution (44) takes on the form 
t/[t + At] = y[t]e- A ^] + l( e - A ^‘] + l).At + 

+ n {e~ A * M (i (tfiM + ei) - h {2C/ 2 [t] + 2e\' 2 (l - ^ /2 )}-At - 

- Af/ 3 [<].At 2 j - ffij-At 2 , (86) 

where 

A¥[<?i] = I (f/iW + ex)- At + i {c/ 2 [t] - 2e\' 2 (l - £i /2 )}'At 2 

is the argument of the exponential. 

Equations (85) or (86) may now be substituted into (84) to give 

£i — y[t + A*] 2 = G[qi], 

and therefore we obtain the implicit relationship 

T[qx] = 0 where T[q x ] = Q[q x ] - q x . 

This equation is easily solved for Qi by the Newton-Raphson method and the solution, when inserted into 
(85) or (86), yields y[t -f At] in terms of y[t]. 

A comparison of the exact solution (79) with an Euler backward difference solution, and with the solutions 
given by equations (85) and (86), with and without the quadratic terms, is given in Table 2. In this table, y^x 
denotes the exact solution (79); ysE the backward Euler solution; y^j the linear implicit solution obtained 
by keeping only the first two terms in (85), and setting f/ 2 = 0; yqn the quadratic implicit solution obtained 
by keeping one extra term in (85) — the first term in the infinite series — and retaining f/ 2 ; yqn the quadratic 
implicit solution obtained by keeping two extra terms in (85) — the first two terms in the infinite series; yix 
the linear trapezoidal (Euler-Maclaurin) solution obtained by keeping only the first two terms in (86), and 
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Table 2: Accuracy of Several Integration Algorithms. 



Value of y at Point A j 

Number of 
Integration 
Steps to A 

VEX 

Vbe 

Vli 

yqn 

yqi 2 

Vlt 

Vqem 

1 

.8230 

.6823 

.7597 

.8107 

.8150 

.8488 

.8247 

2 

.8230 

.7459 

.7800 

.8192 

.8196 

.8331 

.8232 

5 

.8230 

.7895 

.8020 

.8226 

.8226 

.8248 

.8230 

10 

.8230 

.8057 

.8118 

.8230 

.8230 

.8235 

.8230 


Value of y at Point B 

Number of 








Integration 
Steps to B 

Vex 

Vbe 

Vli 

Vqi\ 

Vqi 2 

Vur 

Vqem 

1 

.9895 

.8351 

.9393 

.9677 

.9689 

1.2237 

1.0053 

2 

.9895 

.9154 

.9579 

.9803 

.9808 

1.0450 

.9888 

5 

.9895 

.9630 

.9751 

.9871 

.9872 

.9966 

.9895 

10 

.9895 

.9772 

.9821 

.9888 

.9888 

.9912 

.9895 


setting U 2 = 0; and vqem the quadratic Euler-Maclaurin solution given by (86). The values of y at points 
A and B, where t = 1 and t — 2, were found from the various algorithms when the number of integration 
steps to reach points A and B were taken to be 1, 2, 5, and 10. It may be seen from Table 2 that both the 
linear implicit and trapezoidal approximations using only one step to reach point A, viz. y^i = .7597 and 
Vlt = .8488, are superior in accuracy to the backward Euler algorithm, y^E = .6823. Generally speaking, 
the linear asymptotic solution is as accurate as the Euler backward difference solution with twice the number 
of integration steps. When one extra term is retained in the asymptotic implicit algorithm, the solution using 
one step to reach point A, yQn — .8107, differs from the exact solution, yEx — .8230, by about one percent, 
whilst with one time step the Euler backward difference solution, ysE = .6823, is in error by about seventeen 
percent. There is thus some incentive to retain the second order term in the algorithm, because the accuracy 
achieved in one integration step, with this term included, exceeds that of the Euler backward difference 
algorithm with more than ten times as many integration steps. In all cases, the quadratic Euler-Maclaurin 
solution is the most accurate of the methods considered. 

The results of Table 2 are given graphical interpretation in Fig. 2. Here we plot the relative error incurred 
by the various integration methods versus the number of evenly spaced integration steps used to integrate 
out to point A. The steeper the slope of the curve the greater the rate of convergence. The backward 
Euler and linear implicit methods have slopes of about —1 with linear rates of convergence. The linear 
trapezoidal and quadratic implicit methods have slopes of about —2 with quadratic rates of convergence. 
And the quadratic Euler-Maclaurin method has a slope of about —4 with a quartic rate of convergence. This 
illustrates, but does not provide theoretical proof of, the convergence properties of these methods. There 
appears to be a power-law relationship between the number of integration steps and the error. To check 
this, 1000 integration steps were used with the quadratic Euler-Maclaurin method resulting in an error on 
the order of 10~ 14 , which corresponds to a linear extrapolation of the yQEM curve in Fig. 2. Of note is the 
fact that the Euler-Maclaurin methods have rates of convergence that are twice their order, while the Taylor 
methods have rates of convergence that equal their order. 

Point B is in the asymptotic region where y approaches unity. Here the differences between the Euler 
backward difference algorithm and the present formulations decrease. This is due to the fact that the higher 
order derivatives of U vanish as y approaches its asymptotic value. However, it may be seen from Table 2 
that even at the point B where t = 2, there is still a considerable error in the Euler backward difference 
solution as compared to the present solutions. Qualitatively, the trends present in Fig. 2 for point A are the 
same at point B. Quantitatively, there are differences in the errors for any given method between points A 
and B, but these are slight. 
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The differences between the solutions yqji and yqn are not large, and greater accuracy can therefore 
only be attained by the use of a cubic or higher order solution. Including the third term in the infinite series 
will do little to enhance the accuracy of the quadratic solution. 

When a system of N differential equations is solved, the Euler backward difference solution requires 
the inversion of an N x N Jacobian; whereas, the present formulation generally requires the inversion of 
an M x M Jacobian, where M <C N is the number of independent scalars in the integral solution. This 
is demonstrated in the remaining examples. The advantage of the present method thus manifests itself in 
increased accuracy with a smaller Jacobian matrix. 

9.2 Example 2 

Examples 1 and 4 present differential equations that possess asymptotic limit states in which the second 
and higher order derivatives of U and V vanish. There the derivative U\ increases to a limiting value and 
remains constant until a change in the direction of the forcing function V\ is encountered. In such cases 
the asymptotic approximations provide excellent accuracy. In the present example, which is concerned with 
modelling the spread of the AIDS virus in a susceptible fraction of the population, the application of the 
asymptotic Taylor algorithm results in solutions of the differential equations that are too heavily damped. 

In this example, the function U\ increases at first, but then decays towards zero. Although it is always 
positive, and U is therefore always an increasing function of time, the decay in U\ means that the exponential 
in the integrand does not decay rapidly enough as we depart from the upper limit of the nonhomogeneous 
integral for the asymptotic Taylor methods to provide an accurate representation of its value. Under these 
circumstances, the trapezoidal approximation to the integral given in equation (25) is more appropriate, 
because it weights the contribution to the integral from the upper and lower limits equally. 

Hyman and Stanley [22] have presented a number of models for describing the spread of the AIDS virus 
in a susceptible fraction of the population. In their simplest model, where 

S(t) is the number of susceptible individuals, 

I(t) is the number of infected individuals, 

A(t) is the number of AIDS cases, 

H is the death rate of individuals without AIDS, 

6 is the death rate of individuals with AIDS, 

7 is the rate of developing AIDS of infected individuals, 

i is the probability of infection from an intimate contact with an infected individual, 

c is the average number of contacts between partners, 

r is the average number of new partners per year, and 

So is the susceptible population size before the introduction of the AIDS virus, 

the governing equations are defined to be 


where 


S + (n + A) S = ftS 0 , 

(87) 

i + (7 + /i) I = AS, 

(88) 

A + 6 A = 7 1, 

(89) 

. • I 

X — ic r— 

5 + I 

(90) 


In the simulations that follow the constants have the assumed values: /i = .02, i = .04, c = 1, r = 36, j — .1, 
8 = .5, and So = 10,000,000. 

9.2.1 Implicit Asymptotic Approximation 

For the implicit, linear, asymptotic approximation, the appropriate recursive relations are derived from (19) 
in the form 

A[t + At] = A[t] e~ s ^ + — (1 - , (91) 
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and 


I\t + ATI — lit] Ai + + A<]5[< + A/] / c -(rf/i) A<\ 

l[t + At\-me + ( 7 + Ai )(5[/ + A<] + 7[< + A<])l l )’ 

i c rl[t 4- A<] j 


S[t + At] = S[t] exp 




*S[i + A<] 4- I\t A<] 

pS 0 


+ 


(- 


icrI[t + At] N 
*5[t + At] + 7[t + At] > 


(l - exp - ^ 4- 


» c r/[< 4- At] 

5[t 4- At] 4- I[t 4- At] 


H) 


On putting 


0i = /* + 


icr/[t + At] 

S[t + At] + /[t + At] 1 

rewritten in terms of 

A< 4- — h-e-*' A 
01 


and 


7[f. 



(92) 

0i as 


‘) = S[ ei ], 

(93) 

c _ (7+ ^).aA 3 I[ei] 

(94) 


Substitution of (94) and (93) into (92) gives the required implicit equation for £i, 

_ , icrl{g\] 

=• '‘ + Stoll + /b.] - 

which, after solution with the Newton-Raphson method, can be substituted into the recursive relations (91), 
(93), and (94) to advance the solution from t to t + At. 

9.2.2 Implicit Trapezoidal Approximation 

More appropriate recursive relations can be obtained by employing the trapezoidal approximation in (25) to 
evaluate the nonhomogeneous integral. In this case we write 


rt + At 

ft+At T rt+At 

- (P + A[r]) dr 

+ exp - / (t*+ A[r])dr 

Jr = t 

J(=t [ J T=( 


+ At] = S[t]expl- / (/i + A[r]) dr| + / exp | - / (/i + A[r]) dr \ fiS 0 d£ 

~ S[*] e -*(2j. + A[(]+A[ t + AtJ)A< + 1 ^ Sq ^ +e -A(2^+AM+A[ t + A«]) 

On setting 

01 = ^ + -M* + At]), 

the appropriate recursive relations have the form 

S[t + At] = S[t] e~ e> At + \vSo (1 + e"‘' A ‘) At = 5[p,], 

I[t + At] = I[t] e -(7+rt’ A ‘ 4- § (A [t]S[t] e-^ A< + {2( ft - /i) - A[<]}S[*,]) A< = I[g ,], 


and 


+ A<] = A[t] e-* M + i 7 (/[<] A ‘ + I[e i]) - A* = A[e i], 


(95) 

(96) 

(97) 

(98) 


where Aft + At] has been replaced in terms of Q\ through equation (95). 

The required equation for Q\ is then obtained by substituting (96) and (97) into (90) and (95) in the 
form 


To establish an accurate baseline solution, equations (87-90) were integrated over a period of 40 years 
with an explicit Euler forward difference method employing 8,000 integration steps, and we refer to this 
solution as the ‘exact’ solution for brevity. 

A comparison of the ‘exact’ solution with an explicit Euler forward difference method and the linear 
asymptotic implicit relations in (91-93) employing 400 integration steps each is shown in Fig. 3. We refer 
the reader to the report [22] by Hyman and Stanley for a discussion of the solutions to the differential 
equations. It is evident that the explicit forward difference method provides an underdamped response, 
and that the linear implicit asymptotic solution provides an overdamped response. The implicit trapezoidal 
algorithm given by relations (96-98) effectively averages the amount of under- and overdamping, and a 
comparison of this solution employing 400 integration steps with the ‘exact’ solution is shown in Fig. 4. The 
averaging expectation is borne out, and a very accurate solution is achieved. 

When only 40 time steps are used to integrate the solution out to 40 years the effects are more dramatic. 
As shown in Fig. 4, the trapezoidal solution shows some error. But as shown in Fig. 5, the asymptotic 
solution is heavily overdamped and badly in error, while the explicit forward difference solution is heavily 
underdamped and errs in the opposite direction. 

It appears that the implicit asymptotic algorithm does a good job when U\ is a monotonically increasing 
function. However, when U\ can increase and decrease over a given time increment, the implicit trapezoidal 
(or a higher order Euler-Maclaurin) algorithm is the more accurate. 


9.2.3 Improved Trapezoidal Approximation 

In the preceding approximation, the forcing function on the right hand side of equation (88) is A S. On 
examining equation (90) we find that the forcing function is a function of /, the variable in the differential 
equation. The differential equation can therefore be brought to the homogeneous form, 


i + (7 + ^ 1 = 0, 


and the appropriate trapezoidal recursion relations become 

S[t + At] = S[t] e" + 

+ |/z5 0 (l + e -i( 2 M+A[<]+A[t+Ai])-A«^ At 

and 


I[t + At] = I[t] exp 
On defining the variable Qi as 




27 + 2/i - 


A[t]S[t] A[t + At]S[f + At] 


£1 


"K 


27 + 2 /j 


I[t] I[t + At] 

A[(]5[<] A [f + At]S[t + At] 


H 


m 


/[< + At] 


)• 


we can immediately write 

/[t + At] = J[t]e"*' A ‘ = /[ ei ], 

We can use the definition of A in (90), viz. 

i crl[t 4- At] i c r/[^i ] 


A[t 4- At] — 


S[t -h At] I[t 4" At] «S[t 4" At] 4* 
to eliminate the expression A[t 4- At]S[t 4- A t]/I[t 4- At] from equation (100), and we find that 

A[t]S[t] 


A[f 4- At] = i c v — 27 — 4- 


/[t] 


+ 2^i = A[gi]. 


This expression may then be substituted into (99) to give 

S[t + At] = S[t] e -i(^+M*]+Me,]) A< + I Ai5o ^ + e -i(2M+A[«]+AU,]) _ Af s 


(99) 


(100) 

( 101 ) 

( 102 ) 

(103) 

(104) 
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The recursion relation for A may then be written as 


A[t + At] = A[t] e~ 6 £kt + i jl[t] (e~ s A ‘ + e"« l A ‘) At = Xfo]. (105) 


The required equation for q\ can now be obtained by substituting the expressions for \[qi], I[g i], and S[g\] 
into (100) to give 


^■[ei] = 7+^-| 


/ A[t]5[t] A [ Bl ]S[ei] \ 

l m I[Q 1 ] ) 


- Qi = 0 . 


The solution obtained with the homogeneous version of the equation for I is shown in Fig. 6 using only 
40 time steps out to 40 years. It is clear from this example that much more accurate solutions are achieved if 
the differential equations can be arranged so that the forcing function on the right hand side of each equation 
is independent of the variable being solved for. Thus, we should like the right hand side of the 5 equation to 
be independent of 5, that for I independent of 7, and that for A independent of A. A zero (homogeneous), 
constant, or slowly varying function of the other variables is the ideal forcing function. This is, perhaps, to 
be expected because then the right hand side is known, and we therefore approximate the exact Bernoulli 
conditions of constant U\ and V\ more closely, with all of the nonlinearity then residing in the exponential. 


9.3 Example 3 

The predator-prey equations of Volterra and Lotka [23] form a simple nonlinear system of differential equa- 
tions governing the population of two biological species comprised of predators and prey. These equations 
are known to exhibit closed oscillations similar to the stable vibrations of an undamped pendulum. 

If x denotes the prey population and y the predator population, we can write the coupled system of 
differential equations — in which the constants are taken as unity — in the form 

x — x - xy, (106) 

and 

y-xy-y. (107) 

9.3.1 An Inappropriate Asymptotic Algorithm 

When the preceding equations are rewritten as x + (y — l)x = 0 and y -f (1 — x)y = 0 it is evident that 
the solutions will contain growing and decaying exponential type solutions according as y is less than or 
greater than one for x, and according as x is greater than or less than one for y. It is clear at the outset 
that asymptotic methods are not appropriate here unless we distinguish the ranges of x and y and apply the 
appropriate asymptotic expansion in each region for each equation. We have not done this, but we do show 
that the linear, implicit, asymptotic algorithm — which works well for the first and fourth examples — predicts 
enough damping that the solution (known to perform almost circular stable paths in the x~y plane) spirals 
in towards the steady state solution, x — y — 1, where the populations of predator and prey become fixed 
in size. In their book [23], Thompson and Stewart remark that “One acknowledged deficiency of the Lotka- 
Volterra equations is thus that they have the structural instability of the undamped conservative mechanical 
system, the phase trajectories of which can be topologically changed by the introduction of infinitesimal 
viscous damping.” The implicit asymptotic method, which was shown to exhibit too much damping in 
Example 2, precludes the formation of a stable closed loop in the x~y plane. 

The linear asymptotic recursion relations obtained from the differential equations 



x + yx — x 

(108) 

and 


y+ ly = xy 

(109) 

can be written as 


x[t + Ai] = *[<] e _fl A< + + (1 e"*‘ A ‘) 

(110) 
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and 


y[< + At] = y[t] e At + Q\ x\t + At] (l - e A< ) , 


( 111 ) 


whprp 

ei = y[t + At]. (112) 

Equation (110) can be rearranged to give an explicit dependence on ffi in the form 

x[t + At] = ^ = x[ei], (113) 

1 (1 -e-<> A ‘) 

Q\ 

which may then be substituted into (111) to give 

y[t + At] = y[f] e- A< + (i _ e -**) = (114) 

1 (1 - 

Q\ 

The preceding equation can then be substituted into (112) to yield the required equation for qi 9 viz . 

?[Qi] ~ y[Qi] — ^ 1 = 0 . 

On evaluating the solution by Newton-Raphson iteration, and substituting into (113) and (114), we advance 
from time t to t + At. 

The solution in the x-y plane for this system of equations is shown in Fig. 7 under the initial starting 
conditions x[0] = 1.0 and y[0] = 0.2 using a time step of At = 0.1. It is evident that the solution is 
spiralling in towards the steady state solution x = y = 1 due to the damping that is inherent in the linear 
asymptotic expansion. This spiralling in towards the steady state solution still occurs when much smaller 
time steps are used, but the rate is decreased. If an explicit, Euler, forward difference method or the linear, 
explicit, recursion relation of equation (21) is used to integrate these differential equations with a time step of 
At = 0.1, the solution is underdamped and spirals out to infinity. Figure 8 shows the behavior of the linear, 
explicit, recursion algorithm* With this size time step, a simple Heun — or improved Euler — method [24] 
provides a stable oscillatory response, whilst a fourth order Runge-Kutta method can tolerate even larger 
time steps. In the next subsection we attempt to improve on the implicit asymptotic method, which we 
know to be inappropriate for this problem, by using the implicit trapezoidal approximation which worked 
well in Example 2. 


9.3.2 An Appropriate Trapezoidal Algorithm 

The implicit, trapezoidal, recursion relation in equation (25) can be used to generate the following recursion 
relations from the differential relations in (108) and (109): 

x[t + At] = X [t] e~ i(»W+»[*+ A «])' + I (*[<] «- JM*I+*[*+A*D-A« + x [t + At]) At, 


and 


If we put 


y[t -f At] = y[t] e Ai + \ e + x[t + At]y[t + A t]^ • At. 


e\ = i (yVl + y[t + A *D . 

the recursion relations can be rearranged and expressed in terms of q\ to yield 


x[* + At] = 


.[!]«— ■**(!+ jAl) j , 

i - |a t 


and 


5i!]e-"(l + Hf] Af) _ , , 

*•' + A!l = 


(115) 

(116) 
(117) 
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The governing equation for g x is then found from (115) and (117) as 

?lei] = 5 (»[<] + y[ei]) - ei = 0. 

The variation of the solution in the x-y plane for a time step of At = 0.1 is shown in Fig. 9, where it may 
be seen that the solution does spiral inwards, but at a much smaller rate than the asymptotic solution of Fig. 
7. When smaller time steps are employed the spiralling inwards can be arrested, and a steady oscillatory 
response is obtained. 

The preceding trapezoidal solution was given in order to emphasize the importance of rearranging the 
original differential equations. There is, in general, no unique way of arranging the equations for solution, 
and the ingenuity of the analyst is usually required to obtain an optimum form of the equations for solution. 
In the preceding solution the differential equations were arranged in such a form that X U X — y and 2 U\ — 1 
with forcing functions X V X = x and 2 Vi — xy. As discussed in Example 2, they may also be rearranged into a 
more appropriate homogeneous form with zero valued forcing functions. We obtain the implicit trapezoidal 
solution for this case in the next subsection. 

9.3.3 A Better Trapezoidal Approximation 

It is possible to rearrange the differential equations (108) and (109) into the homogeneous form 

x + (y - \)x = 0, 

and 

y + (1 - x)y = 0. 

If we now let 

ei — \ (y[<] + y[< + At] - 2) and £2 = \ (2 - a:[f] - x[t + At]) , 

then the implicit, trapezoidal, recursion relations with no forcing terms are obtained in the form 

x[t + At] = x[t] e~ erAt and y[t + At] = y[t] (118) 

The relations for g x and g 2 can then be written as 

8i = i(»[<](l + «-•»•**) -2) and £2 = \ (2-*[t](l + e"'>- A ‘)) , (119) 

which when combined gives the required equation for g x , i.e. 

F[e\] = 5 {y[<] (l + ex p[~| {2 — x[t] (l + e _<?1 A ')} -At]) - 2} - Qi = 0. 

On solution, g x may be substituted into (119) and thence into (118) to step from t to t + At. 

Results from this solution procedure are shown in Fig. 10 for a time step size of At = 0.1, where the 
improvement is at once apparent. This solution still gives a tolerable but ragged stable loop when the time 
step is increased to At = 0.3, a size at which the fourth order Runge-Kutta solution tends to show an even 
more ragged approximation to the stable solution. The improved accuracy achieved when the right hand side 
of each equation can be made homogeneous, constant, or weakly dependent on the other solution variables 
is once again apparent. 

9.4 Example 4 

In the previous examples the Jacobians are of dimension 1, and therefore the implicit algorithms developed 
for their solution are straightforward. In this example we present a viscoplastic model where the governing 
system of differential equations is sufficiently complex to warrant a Jacobian of dimension greater than 1. 
Our implicit integration algorithms have an important advantage over classical algorithms — like backward 
Euler integration — in that they have the distinct possibility of a substantial reduction in the order of the 
Jacobian. In the AIDS model of Example 2, there is a reduction in order from a 3 x 3 to a 1 x 1 for the 
Jacobian matrix; for the predator-prey model of Example 3, there is a reduction in order from a 2 x 2 to 
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a 1 x 1 for the Jacobian matrix; and for the viscoplastic model considered below, there is a reduction in 
order from a 13 x 13 to a 2 x 2 for the Jacobian matrix. This reduction in order of the Jacobian matrix is 
of significant importance from a computational viewpoint. 

Viscoplastic constitutive models [7, 19] are currently finding application in describing the rate dependent 
plastic response of metallic structures that undergo significant temperature change over their duty cycle. We 
developed the viscoplastic model [25] presented below for metallic polycrystalline materials. It is presented 
here in a simplified form that is valid for temperatures in excess of about one half the melting temperature. 
In this model, 
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is 
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T 


is 

the 

Y 


is 

the 
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is 
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E%j — €%j - 

- 3 £kk$ij 
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the 



is 
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is 
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The isotropic, stress/st rain, 


where 


drag strength, 
absolute temperature, 
yield strength, 

(deviatoric) back stress tensor, 
deviator ic strain tensor, 

deviatoric stress tensor, 
infinitesimal strain tensor, 

(deviatoric) plastic strain tensor, and 
Cauchy stress tensor. 

constitutive equations are given by Hooke’s law, 

Sij = 2p (Eij-e ?•), 

<Tkk = 3 k (£** — a AT 6kk) > 


AT = T - To 
a 

K 


% 


is the temperature change with respect to a reference, To, 

is the coefficient of thermal expansion, 

is the bulk modulus, 

is the shear modulus, and 

is the Kronecker delta; 1 if i = j t otherwise 0. 


( 120 ) 

( 121 ) 


Here one sees that the plastic and thermal strains, and ot AT6 |; -, are eigenstrains for the deviatoric and 
hydrostatic responses, respectively. 

The evolution of plastic strain is described by the viscoplastic flow law, 



*-*11*11 iieV 

(122) 

with a kinetic law, 





(123) 

and two evolutionary laws, 

*="(*" Jot "'l 



) > (124) 


Y = v(h[Y] llel-Kr.y]) 

1 , (125) 


where 


h is the nonlinear hardening parameter for yield strength, 

H is the hardening modulus for back stress, 

L is the limiting state for the dynamic recovery of back stress, 

r is the thermal recovery parameter for yield strength, 

Z is the Zener parameter, 

t i is the hardening modulus for yield strength, 
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? 


= Sij — Bij 

¥ p \\ = y /^A 

||E||= -^EyEy 
(x) 


is the thermal diffusivity, 
is the effective stress, 

is the magnitude of plastic strain rate, 
is the magnitude of effective stress, and 

is the Macauley bracket of x whose value is x whenever x > 0, otherwise it is 0. 


Viscoplasticity is an internal state variable theory where the internal state variables, B ij and Y in this 
case, are governed by separate evolution equations. The evolution of internal state follows a competitive 
process between hardening and recovery (both thermal and dynamic) mechanisms. The back stress, Bij, is a 
phenomenological representation of an internal stress state; it reflects the stress fields set up by dislocations 
in their heterogeneous substructures. The yield strength, Y , is a phenomenological representation of material 
strength; it reflects the density of dislocations. 

We will consider a viscoplastic model characterized by material functions that are power-laws: 


m = 

r q l 
exp r«rJ’ 

f(lis||-y>] 


f<nsii-m 3 

D 


V D ) ’ 

L[Y] = 

(H— (£f 



f Y \ 3 -" 

h[Y] = 

(»c) ’ 

r[T,Y] = 

'«(«'■ 


where 

C is the creep strength, 

n is the creep exponent, 

Q is the activation energy for creep (or self diffusion), 

R is the universal gas constant, 8.314 J/mol-K, and 
y is the fraction of yield strength to applied stress at steady state. 


Because of the chosen forms for the material functions L, h , and r, this viscoplastic model analytically 
reduces to the classical theory of creep [26] at steady state, where there is no evolution of the internal state 
variables, and where 


= km 


\ C ) ||S|| 


defines the creep rate, with ||S|| = y \SijS{j characterizing the magnitude of deviatoric stress. 
For illustrative purposes, we will consider the following set of elastic material constants: 


a = 20 x 10~ 6 °CT 1 , 

k = 83, 750 MPa, 

/i = 30,000 MPa, 


and viscoplastic material constants: 


C = 0.8 

MPa, 

D = 0.016 

MPa, 

H=n/2 

MPa, 

n = 5, 


Q = 200,000 

J/mol. 

y =0.1, 


T] = fi 

MPa, 
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which characterize copper in the neighborhood of 500°C. 

If we differentiate equation (120), combine it with (122), and subtract (124) from both sides; and then 
rearrange equations (124) and (125), the following set of differential equations is obtained: 


• fo+£)J!£]i _ 2 A ... H He'll 

+ iTFTi ^ Fr\7i > 


Y + rj 


E|| ■ L[Y ] 

iig p ir 


^[r, 


-) Y = 0, 


■ ^lie'll _ ^H^ll 
,j+ L[Y] ij ~ || E || 


This system of equations has the form 

Xa + X a = & 


(for a = 1, 2, 3), 


(126) 

(127) 

(128) 


where 


and 




x, = y. 


X 3 = Btj, 


(fi + H)\\e»\ 

lUl ~ — imt - 


iU x 




r[r,y] -ft[y]||e P || y 


3^1 


#l|e p ll 

i[y] 


iKi = 2fiEij + 


£J1£J1 b 

i[y] 


2 V) = 0 , 


3 K 1 


IIS || 


It is now apparent why the Jacobian is at most a 3 x 3 matrix (one for each of the three differential equations) 
whereas the Jacobian for backward Euler integration is a 13 x 13 matrix (one for each of the thirteen spatial 
dimensions: six for the effective stress, six for the back stress, and one for the yield strength). However, 
the Jacobian must be further reduced to a 2 x 2 matrix. This occurs because \U\ and 3 C /1 both become 
asymptotically proportional to ||£: p || at steady state, and are therefore linearly dependent at steady state, 
Hence, the evolution equation for back stress (128) must not contribute to the construction of the Jacobian 
in order to prevent the Jacobian from becoming singular. 


9.4.1 An Asymptotic Algorithm 

For linear, implicit, asymptotic integration, one has 


X Q [« + A<] = X a [<]e-^ At + 


aVl [Xy[t + A*] 


Qa 


(1-e- 


Qc Af\ 


for the integration algorithm, where 

(A<+tf)||£ p | 


ei 


and Q 2 








which are both evaluated at time t + At. Thus, the set of iteration functions becomes 


-^1 [{^}a] = 


0 i + H) || g P [{gl, g2}A]|| 

l|S[{ei>A)|| 


{ei>A, 


and 

^(r[T,y[{^h]]-h[y[{e2}A]l ||e p [{e 1 ,52}A]||) 

y[{e 2 }A] - {^}a, 

which is solved through Newton- Raphson iteration. The derivatives required to construct the Jacobian in 
equation (70), e.g. d\\i p \\/dgi, etc., have been determined numerically, for in this case the determination 
of numerical derivatives is computationally more efficient than determining them analytically. 
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The Newton-Raphson iterations are accomplished by the following scheme. First, values of Q\ and g 2 
are guessed. In addition, we also guess the values of E and Y. The forcing vectors, iV\ } 2^1 > and 
3V1 are then computed and used with q\ and g 2 in the recursion relationships to determine X \ , X 2 , and 
A3. Estimates of E;y , B tJ -, and Y are now available to compute improved estimates of the forcing vectors, 
\V \ , 2V1, and 3V1 for the next iteration. It is evident from the preceding that this approximate algorithm is 
basically one iteration in arrears in estimating E fJ -, , and Y. 

The capability of this integration method is demonstrated in Fig. 11. The solid line was obtained using 
250 integration steps, and may be considered to be the exact solution. In this example, a homogeneous block 
of material is sheared in a specified direction at a constant rate of straining to a fixed value of engineering 
strain, i.e. 7 = 2ei2- This block of material is then sheared in the opposite direction to a self-similar value 
of engineering strain. The direction of loading is changed one last time to complete the loading cycle. The 
increase in the shear stress from the first to the third reversal is due to a gradual increase in the value of the 
yield strength, Y, over the loading history, i.e. the material is gradually getting stronger because of work 
hardening. The smooth curvature observed in the stress-strain response of each of the three loading segments 
is due to the rapid evolution of the back stress, , over each segment. For a more detailed discussion of 
the model, see Freed and Walker [25], and of viscoplasticity in general, see Lemaitre and Chaboche [7]. The 
differences between the exact response and those predicted using 25 and 3 integration steps are very small; 
however, there is a measurable error for the case of 10 integration steps at the knee of the curve. At first 
glance, the fact that 3 integration steps does better than 10, seems contradictory. However, the contradiction 
is apparent only, because the three integration points are close to their asymptotic solutions, where the 
implicit, asymptotic, integration method is very accurate. The regions within which the ten integration 
points are in greatest error are the transient domains where the back stress is rapidly evolving. The use of 
a higher order, asymptotic, integration method would reduce this error. One important observation about 
this integration method is that the error generated in the transient domains does not propagate with the 
solution into the asymptotic domains. This is because the correct asymptotic solution of the differential 
equation is contained within the linear, implicit, asymptotic, integration method. 


9.4.2 A Trapezoidal Algorithm 

For implicit, trapezoidal, Euler-Maclaurin integration, one has 

X a [< + A<] = X a [f]e- e " At + | A * + a v\ [x 7 [f + A<]j) -At 

for the integration algorithm, where 


/||£W|i P p [< + A<]j| \ 

* 2 VllSWH ||E[< + A<]||y ’ 

and 

_ *? fr[t] -h[t] || e p [<] || r[t + At] — h\t + At] ||e p [t + At] || 

— 2 \ y[<] + y[< + A<] 

The set of associated iteration functions is therefore given by 


and 


T \ [{^7 U] = 


fi + H 
2 


f CTH , ||g P [{gl.g2}A]|| \ , , 

viiEwir m{ ei }x}\\ ) ^ ljA 




77 l r[£] - h[t] ||£ p [<]|| 

2 ^ Y[t] 


r 

4- — 

T,Y[{ ei}x] 

- h 

Y[{e 2 }a] 

l|e P [{f?l,(?2h]|| > \ 





{Q2}X, 


which is solved through Newton-Raphson iteration, as described in the previous subsection. 
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The capability of this integration method is shown in Fig. 12 for the same loading history that was used in 
Fig. 11. The solid line represents 250 integration points, and is the exact response. For all practical purposes, 
it is identical to the solid line of the previous figure. Good predictions are made with 50 integration points. 
When 25 integration points are used, there is an error of about six percent between the exact and predicted 
values. As the time steps are increased in size this error increases. For 10 time steps the error approaches 
fifty percent. Clearly one must monitor the size of the time step when using trapezoidal integration so as 
to secure accurate answers. It is also apparent in Fig. 12 that the trapezoidal solution does not converge 
toward the asymptotic solution, as is the case in the linear, asymptotic, integration method, for the reasons 
stated below. 

The inaccuracies attendant on the larger time steps are due to the representation in equation (126). The 
right hand side of this equation contains the term \\i p \\ which depends in a highly nonlinear manner on the 
variable E |; * on the left hand side of the equation. As previously demonstrated in Examples 2 and 3, better 
accuracy can be achieved if this circumstance is avoided. 

10 Concluding Remarks 

Generally, the methods which employ asymptotic Taylor expansions to evaluate the nonhomogeneous integral 
are suitable for stiff and boundary layer problems, whilst the trapezoidal and asymptotic Euler-Maclaurin 
expansions of the integral are more suited to general equations. The integration algorithms developed herein 
offer marked improvements in accuracy and stability over existing integration algorithms for systems of 
nonhomogeneous, nonlinear, first order, ordinary differential equations. These new algorithms are reviewed 
in the Overview, section 3. Several of these algorithms are asymptotically correct, thereby enabling large 
time steps to be used while retaining stability and accuracy. In all cases, accuracy is greatly enhanced if 
the right hand side of the equations (the nonhomogeneous contributions) can be made to be zero valued 
(homogeneous), constant valued, or slowly varying functions of the other variables in the system of equations 
being integrated. Four examples are presented to demonstrate the viability of these new algorithms. Various 
solution strategies are given for each example to demonstrate the advantage of proper construction of a 
solution procedure. 

11 Acknowledgment 

This work (for K.P.W.) was supported by the NASA-Lewis Research Center under Grant NAG3-882, and 
by the Department of Energy under Contract DE-AC02-88ER13895. 


References 

[1] C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations , (Prentice-Hall, 
Englewood Cliffs, NJ, 1969). 

[2] J.-W. Ju, /. Engr . Mech, 116, 1764 (1990). 

[3] K. P. Walker, in Advances in Inelastic Analysis , AMD-Vol. 88 & PED-Vol. 28, edited by S. Nakazawa 
et a/., (ASME, New York, 1987), p. 13. 

[4] S. C. Dennis, Proc. Camb. Phil. Soc 56, 240 (1959). 

[5] A. D. Freed and K. P. Walker, in Visco-Plastic Behavior of New Materials, PVP-Vol. 184 & 
MD-Vol 17, edited by D. Hui, and T. J. Kozik, (ASME, New York, 1989), p. 1. 

[6] A. Chulya and K. P. Walker, A New Uniformly Valid Asymptotic Integration Algorithm for Elasto- 
Plastic-Creep and Unified Viscoplastic Theories Including Continuum Damage , (NASA-Lewis Research 
Center, TM-102344, 1989). 

[7] J. LEMAITRE AND J.-L. CHABOCHE, Mechanics of Solid Materials , (Cambridge University Press, Lon- 
don, 1990), ch. 6. 


36 



[8] E. L. INCE, Ordinary Differential Equations , (Dover, New York, 1956), pp. 21 and 531. 

[9] J. BERNOULLI, Acta Eruditorum publicata Lipsce , 113 (1697). 

[10] R. W. Hamming, Numerical Methods for Scientists and Engineers, 2 nd edition , (Dover, New York, 

1986) , p. 421. 

[11] P. S. Laplace, Theorie Analytique des Probability , Vol 7, (Paris, 1820). 

[12] N. Bleistein and R. A. Handelsman, Asymptotic Expansion of Integrals, (Dover, New York, 1986), 
ch. 5. 

[13] E.T. Copson, Asymptotic Expansions , (Cambridge University Press, London, 1965). 

[14] A. Erd^LYI, Asymptotic Expansions , (Dover, New York, 1956). 

[15] N. G. DE Bruijn, Asymptotic Methods in Analysis , (Dover, New York, 1961), pp. 40-41. 

[16] H. B. Dwight, Tables of Integrals and Other Mathematical Data, 4 th edition , (Macmillan, New York, 
1961), p. 134. 

[17] K. Knopp, Infinite Sequences and Series , (Dover, New York, 1956), pp. 111-124. 

[18] J. W. L. Glaisher, Quart J . Math., 14, 79 (1875). 

[19] — , Unified Constitutive Equations for Creep and Plasticity, edited by A. K. Miller, (Elsevier, New York, 

1987) . 

[20] T. Y. P. Chang and I. Iskovitz, private communication, (University of Akron, 1990). 

[21] R. D. Krieg, in Transactions of 4 th International Conference on Structural Mechanics in Reactor 
Technology , (SMIRT IV, San Francisco, 1977), paper M6/4. 

[22] J. M. Hyman and E. A. Stanley, Using Mathematical Models to Understand the AIDS Epidemic, 
(Los Alamos National Laboratory, Report LA-UR-87-3078, 1987). 

[23] J. M. T. Thompson and H. B. Stewart, Nonlinear Dynamics and Chaos , (John Wiley and Sons, 
New York, 1988), pp. 44-47. 

[24] W. E. Boyce and R. C. DiPrima, Elementary Differential Equations and Boundary Value Problems, 
2 nd edition , (John Wiley and Sons, New York, 1969), p. 346. 

[25] A. D. Freed and K. P. Walker, A Viscoplastic Model With Application to LiF-22%CaF 2 Hypereu- 
tectic Salt , (NASA-Lewis Research Center, TM 103181, 1990). 

[26] F. K. G. Odqvist, Mathematical Theory of Creep and Creep Rupture, 2 nd edition, (Oxford University 
Press, London, 1974), ch. 5. 


37 



t 


Figure 1: Exact solution and points of integration for Example 1. 



Number of Integration Points 

Figure 2: Relative error of integration compared with the number of steps used to integrate out to point A 

in Example 1. 
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Explicit, forward Euler and linear, implicit, asymptotic (Taylor) solutions for the AIDS model of 
Example 2. 400 integration points. 




Figure 4: Implicit trapezoidal (Euler-Maclaurin) solutions for the AIDS model of Example 2. 
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Figure 6: 


Homogeneous, implicit, trapezoidal (Euler-Maclaurin) solution for the AIDS model of Example 2. 
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Engineering Shear Strain 

Figure 11: Linear, implicit, asymptotic (Taylor) solution for the viscoplastic model of Example 4. 7 = 
±0.001 s- 1 . T = 500°C. Y| (=0 =l MPa. 
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Figure 12: Implicit trapezoidal (Euler-Maclaurin) solution for the viscoplastic model of Example 4. 7 = 
±0.001 s" 1 . T = 500°C. Y| t=0 = 1 MPa. 
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